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I  I  I 


ABSTRACT 

The  modeling  of  a  physical  phenomenon  is  basic  both  to 
science  and  engineering.  While  little  scientific  thought  has  gone 
towards  the  problem  of  model  formulation,  substantial  effort  has  been 
devoted  to  the  process  of  selecting  the  best  model  from  among  the 
competing  models  to  describe  a  physical  phenomenon.  It  is  intended 
here  to  introduce  some  improvements  and  standardize,  the  statistical 
model  discrimination  procedure. 

Present  expressions  for  the  expected  likelihood  of  data 
involve  working  with  the  entire  set  of  data.  New  expressions  are 
developed  that  utilize  only  the  sufficient  statistics  of  the  data 
thereby  affecting  substantial  savings  in  computation  and  avoiding 
large  round  off  errors.  Justification  of  the  conventional  least 
squares  follows  as  a  special  case  of  the  Bayesian  method  from  these 
new  developments. 

Frequently,  it  is  necessary  to  perform  further  experiments 
to  discriminate  between  several  rival  models.  Hitherto,  no  exact 
criterion  has  been  available  for  the  sequential  design  of  experiments 
so  as  to  obtain  maximum  information  in  picking  the  right  model.  An 
expected  entropy  change  criterion  in  line  with  Shannon's  concept  of 
entropy  as  a  function  of  information  about  the  state  of  the  system, 
is  presented  here. 

Finally,  certain  practices  of  Bayesian  procedure  are  criti¬ 
cally  examined  from  the  conceptual  viewpoint.  Examples  of  model 
building  in  Petroleum  Engineering  and  kinetics  are  presented  in  the 
last  chapter. 
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1.  INTRODUCTION 


1 . 1  Mathematical  Models  for  the  Real  World 

A  mathematical  model  is  an  abstract  formulation  aimed  to 
simulate  the  mechanism  of  a  real  world  phenomenon.  A  true  mathe¬ 
matical  model  accounting  completely  for  the  physical  phenomenon  is 
a  hypothetical  concept  and  is  seldom  realized  in  non-trivial  situ¬ 
ations.  Nevertheless,  it  is  a  useful  tool  for  the  scientists  and 
engineers  in  representing  the  important  characteristics  of  the 
actual  mechanism.  An  analyst  is  often  concerned  with  a  situation 
where  a  dependent  variable  y,  is  related  to  a  set  of  independent 
variables  x_  =  (x^ ,  ,  ...»  x^).  For  example,  a  stock  broker 

might  be  studying  the  dependency  of  a  company's  share  value  y,  on 
the  total  investment  of  the  company  x^ ,  the  capital  investment  rate 
per  year  x^ ,  the  company's  total  annual  sales  x^ ,  the  bank  lending 
interest  rate  x^,  etc.  Theoretically,  there  exists  a  relation  of 
the  type 

y  =  <J)  (x ,  a_)  (1.1) 

between  the  dependent  variable  y  and  the  set  of  independent  variables 

x,  where  a  =  (a.,  aOJ  ...,  a  )  is  the  set  of  m  parameters  which  may 
—  —12m 

be  the  physical  constants  of  the  system  under  study.  Usually,  a 
simpler  relation  that  closely  approximates  the  outstanding  features 
of  (1.1)  is  satisfactory. 

A  mathematical  model  that  describes  the  mechanism  reasonabl 
well  might  also  be  useful  to  the  scientists  in  compressing  a  large 
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amount  of  experimental  data  into  a  single  structure.  In  addition, 
the  model  may  provide  new  insight  to  the  phenomenon,  suggest  further 
experiments  and  improvements  in  the  model,  the  purpose  here  being 
to  establish  the  true  nature  of  the  process.  For  an  engineer,  on 
the  other  hand,  the  quantitative  description  is  usually  required 
for  predictive  purposes.  The  mathematical  model  is  a  tool,  where¬ 
by  the  engineer  can  employ  analytical  techniques  in  the  solution 
of  his  problems.  With  the  help  of  the  model,  he  can  isolate  the 
important  aspects  of  the  problem  and  investigate  alternate  designs. 

However,  all  models  do  not  have  to  be  mechanistic  models. 
Depending  upon  the  purpose  of  the  investigation,  an  empirical 
relation  might  be  entirely  satisfactory.  It  would  be  a  waste  of 
resources  to  try  to  obtain  a  mechanistic  model  where  the  response 
to  the  set  of  independent  variables  is  quite  smooth  and  the  only 
purpose  is  to  have  a  working  design  in  the  region  of  interest  ( R I ) , 
particularly  when  extensive  data  is  available  in  tfie  region.  A 
polynomial  type  empirical  model  of  the  form 

m 

y  =  l  a.  .cf>.  (x. )  (1  .2) 

i  =  l 

will  generally  be  satisfactory  and  all  that  is  required  is  to  furnish 
a  suitable  routine  to  generate  successive  functions  to  achieve  the 
desired  accuracy.  Response  surface  studies  deal  primarily  with  the 
empirical  formulation  of  the  models  [see  for  example,  Davies  (1953), 
Box  (1954,  I960,  1964),  Hunter  (1958,  1959a,  1959b)]. 
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In  other  situations  where  the  primary  interest  is  not 
merely  to  predict  the  response  over  a  limited  region  but  rather 
to  elucidate  the  mechanism  or  to  obtain  meaningful  results  in 
regions  relatively  unknown,  the  concern  of  the  investigator  is 
to  obtain  a  mechanistic  representation.  This  kind  of  model  can 
provide  useful  information  where  the  optimization  of  design  leans 
heavily  on  the  knowledge  of  what  actually  is  taking  place.  The 
mechanistic  models  themselves  are  seldom  exact,  and  the  results 
in  far  off  regions  may  not  be  reliable.  However,  an  empirical 
model  lends  absolutely  no  assurance  to  extrapolation. 

The  subject  matter  of  this  thesis  is  concerned  mainly 
with  the  technique  of  obtaining  the  best  possible  mechanistic 
model.  The  objectives  are  to  devise  methods  for:  (a)  obtaining 
the  posterior  model  probabilities,  (b)  design  of  experiments  for 
discriminating  among  rival  models,  (c)  examination  of  the  consis¬ 
tency  of  the  existing  Bayesian  model  discrimination  techniques. 

1 .2  Model i ng  Process 

The  utility  of  the  mathematical  models  to  the  scientists 
and  engineers  make  the  modeling  process  itself  an  important  operation. 
Smallwood  (1968)  illustrates  the  conventional  iterative  process  of 
model  building  as  shown  in  Figure  1(a).  The  first  step  is  to  formu¬ 
late  the  requirements  of  the  model.  As  we  have  seen,  if  the  interest 
of  the  engineer  is  merely  to  have  an  approximate  response  over  a 
finite  region  of  interest,  a  simple  empirical  model  will  do.  Alter¬ 
natively,  resort  to  a  complex  mechanistic  model  might  be  necessary. 


FIGURE  1.1(a)  MODELING  PROCESS 
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A  trade-off  between  the  computational  complexity  and  the  accuracy  of 
the  representation  is  essential  during  the  course  of  this  step. 


The  second  step  involves  hypothesizing  the  form  of  the 
model.  This  is  the  most  crucial  step  since  the  form  of  the  model  is 
directly  related  to  the  physical  mechanism.  A  close  consultation 
between  the  statistician  and  the  scientist  or  engineer  is  required. 

The  difficulty  might  be  overcome  by  training  the  engineer  to  be 
his  own  statistician.  Once  a  form  is  ascribed  to  the  model,  the 
next  step  is  to  estimate  its  parameters  against  the  available  data. 

The  model,  then  is  tested  against  the  original  requirements.  An 
alternate  model  form  can  be  tried  if  the  requirements  are  not  met. 

Deciding  on  the  form  of  the  model  as  in  step  two  above, 
is  to  recognize  the  salient  features  of  the  process  and  ?s  called 
process  identification.  The  subsequent  step  of  ascribing  certain 
values  to  the  parameters  of  the  model  is  termed  parameter  estimation. 
The  process  of  selecting  the  best  parameter  values  can  be  viewed 
more  concretely  by  developing  a  cost  function  of  the  parameter  values. 
The  objective  function,  then  is  to  minimize  the  cost  of  using  a 
particular  set  of  parameters  in  the  design. 

Given  a  model  form  (pj  is  decided,  the  total  cost  of  using 
a  particular  estimate  a_  for  the  set  of  parameters  can  be  expressed  as 


Z  (a. )  =  f  (a. )  .Z  (a. ,  a. ) .  da_. 

J  J  J  J  J  J 


(1.3) 


where  f  (a_. )  is  the  joint  distribution  of  the  parameters  a_.  given  the 
J  J 

model  form  ^  is  true  and  Z  (£_.  ,  a_)  is  the  cost  of  using  a_  as  an 
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estimate  of  the  parameters  when  the  actual  value  is  a. 

The  cost  of  using  a  particular  set  of  parameters  may  be 
interpreted  as  the  savings  that  could  be  realized  in  terms  of 
eliminating  overdesign  had  the  parameter  values  been  known  for 
certain.  This,  then  represents  the  additional  amount  that  could  be 
spent  for  having  the  complete  information  about  the  system.  Stated 
differently,  it  is  the  price  one  would  be  prepared  to  pay  for 
eliminating  uncertainty  associated  with  the  situation  through  addi¬ 
tional  data  acquisition.  The  economic  balance  between  the  cost  of 
the  experiments  and  the  savings  realized  in  the  design  due  to  reduced 
uncertainty  is  the  key  in  deciding  the  extent  of  experimentation. 

However,  there  frequently  arise  situations,  where  many 
model  forms  describe  the  phenomenon  about  equally  well.  In  such 
cases,  it  might  be  desirable  to  consider  all  such  model  forms  simul¬ 
taneously,  with  each  model  form  having  a  set  of  unknown  parameters. 

This  is  specifying  a  model -space  called  a  metamodel  [Smallwood  (1968)] 
which  contains  all  eligible  models,  each  with  all  possible  values  of 
the  parameters.  Process  identification  and  parameter  estimation 
are  now  handled  as  a  joint  problem.  The  termination  of  the  model 

discrimination  process,  in  this  case,  is  decided  by  the  cost  function, 

that  is  the  expectation  of  the  cost  function  associated  with  all  the 
m  individual  model  forms. 


Z 


m 


l  P':  J  f  (l; )  <  L(y  |a_  )  .Z(£  ,  a_)  .da_ 

j  =  l  j  -00  J  J  J  J 


m 


l  P't  I  f  (a.).L(^|a_.).da_ 
j  =  ]  J  _oo  J  J  J 


(1.4) 
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where  L  ( y_  |  a_^ )  is  the  likelihood  of  observing  the  data  given  the  j  th 
model  form  is  true  with  parameters  a_.  and  pV  is  the  prior  probability 
of  the  jth  model  being  true.  Figure  1.1(b)  illustrates  the  modeling 
process  with  the  metamodel  approach. 

The  concept  of  a  metamodel  has  been  used  in  Chapter  3  to 
develop  a  simple  and  useful  expression  for  the  discriminating  function. 

1 . 3  Review  of  the  Previous  Work  in  Model  Discriminations 

The  model  discrimination  procedure  regarded  as  the  process 
of  selecting  the  best  model  from  the  various  eligible  models,  is  well- 
established.  Considerable  work  has  been  done  in  this  field.  After 
the  preliminary  screening  of  the  data  to  eliminate  the  poor  models, 
model  discrimination  may  be  categorized  in  the  following  three  areas: 
(a)  parameter  estimation,  (b)  computing  the  model  discrimination  index- 
likelihoods  of  the  data,  and  (c)  sequential  design  of  experiments. 

These  areas  will  be  reviewed  briefly  to  serve  as  an  introduction  to 
the  text  of  this  thesis, 

1.3-1  Review  of  the  Previous  Work  in  Model  Discriminations 

A  comprehensive  treatment  of  parameter  estimation  of  linear 
and  non-linear  models  is  given  by  Lindley  (1965),  Graybill  (1961)  and 
Raiffa  and  Schlaifer  (1961).  Lee  ( 1 968 )  and  Donnelly  and  Quon  (1970) 
have  described  parameter  estimation  procedures  when  the  models  are  in 
the  form  of  differential  equations. 

The  Bayesian  approach  to  parameter  estimation  has  direct 
advantages  over  the  least  squares  method  in  that  the  prior  knowledge 
about  the  parameters  can  be  effectively  utilized  in  the  Bayesian  method 
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FIGURE  1.1  (b)  METAMODELING  PROCESS 
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For  example,  the  beliefs  about  the  parameters  of  a  model  may  be 
expressed  in  the  form  of  prior  estimates  of  the  parameters  and 
a  prior  variance-covariance  matrix  C^.  The  prior  distribution  of 
a  may  then  be  written  as 


f 1  (a)  =  f  (a  |u  ,  C  ) 

“  N  —  —O  — O 


=  (2tt)  ?/2.  |  C  |  !/2.  exp{-  [a-u  ]lC  1  [a-y  ] } 


0.5) 


Now  if  the  errors  are  assumed  to  be  normally  distributed  with 

2 

mean  zero  and  known  variance  a  ,  the  posterior  distribution  of  the 
parameters  after  n  observations  may  be  approximated  by  the  multivariate 
normal 


where 


and 


f"(ii =  VaJv  V 


Hn 


“O 


C_,y 
l-o  -o 


+ 


-1 


(1.6) 


X_  being  the  n  x  m  design  matrix  of  the  partial  derivatives  of  (x_,  a) 
with  respect  to  the  'm*  parameters  a_ *  (a^,  a^,  ...»  a^) .  An  element 
of  X  i s  denoted  by 


[  (xj  *a) / 9^a j ] ^  ^ 

—  — o 


i  =  1 ,  2,  .  .  .  ,  n 

j  ~  1»  2 ,  •  •  •  i  ^ 


(1.7) 


where  a  is  a  guessed  value  of  a  and  x.  is  the  experimental  condition 

- O  —  ”1 

for  the  i th  observation. 
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In  addition  to  i ncorporat i ng  the  prior  knowledge  about  the 
parameters,  the  Bayesian  method  gives  direct  reliability  estimates 
of  the  parameters  in  terms  of  the  posterior  variance-covariance  matrix. 

For  the  case  when  error  variance  is  not  known,  a  joint 
estimation  of  the  variance  and  the  parameters  can  be  obtained  by 
Bayesian  analysis  [see  Raiffa  and  Schlaifer  (1961),  Quon  (1970)]. 
However  some  sort  of  estimate  is  obtained  for  the  error  variance  and 
in  subsequent  calculations,  this  estimate  is  treated  as  if  it  were 
a  known  quantity.  The  relative  utility  of  the  various  estimates  of 
the  error  variance  is  examined  in  Chapter  4,  from  the  view  point  of 
model  discrimination. 

Sometimes  the  investigator  may  want  to  determine  the  param¬ 
eters  of  a  model  with  the  maximum  possible  precision.  Box  and  Lucas 
(1959)  have  proposed  a  criterion  to  select  n  data  point  settings  that 
will  give  the  most  precise  parameters.  The  criterion  maximizes  the 
determinant  of  X^'X.  The  maximization  of  the  determinant  of  X^X. 
corresponds  to  the  minimization  of  m-d i mens i ona 1  volume  of  the 
confidence  region  in  the  parameter  space.  Directly,  it  leads  to  a 
minimal  contribution  in  the  variance  covariance  matrix  of  the  model, 
thereby  increasing  the  precision  of  the  parameters.  A  chemical  engin¬ 
eering  problem  has  been  studied  using  the  minimum  volume  design 
principle  by  Kittrel,  Hunter  and  Watson  (1966). 

A  Bayesian  procedure  for  estimating  common  parameters  from 
several  responses  was  developed  by  Box  and  Draper  (1965)*  The  method 
is  shown  to  give  more  precise  estimates  of  the  parameters  than  the 
conventional  least  squares.  Draper  and  Hunter  (1968)  developed  a 


sequential  design  criterion  to  obtain  the  most  precise  estimate  of 
parameters  for  more  than  one  response. 


1.3.2  Model  Likelihoods 


Least  sum  of  error  squares  is  the  discriminating  criterion 
for  the  common  least  squares  technique.  A  model  form  that  gives  the 
lowest  least  sum  of  error  squares  is  taken  to  be  the  true  model  and 
the  remaining  models  are  discarded.  Justification  of  favouring  a 
particular  model  by  this  primitive  method  is  given  in  Chapter  2 
from  a  Bayesian  statistical  analysis. 

it  seems  more  logical  to  compute  the  probability  (called 
the  likelihood)  of  the  occurrence  of  the  data  given  that  a  particular 
model  is  true,  and  then  draw  inferences  from  these  likelihoods.  The 
likelihood  of  observing  the  given  data  conditional  upon  jth  model 
being  true  with  parameters  a^  is  given  by  [Raiffaand  Schlaifer 
(1961),  Lind  ley  (1965)] 


L(^Ja.,  $.)  =  (  2ttc  : 2 )  n^2-  exp.{-  ^-v} 

J  J  J  2a.2 


1  1-  I-  tL- 

exp{  -  —  [a.  -a.  ]  [a. -a_  ]  }  (1.8) 

J  J  a  #  J  J 

J 

2 

where  n  is  the  number  of  observations,  a.  is  the  error  variance  for 

J 

the  jth  model  and  includes  both  the  measurement  error  and  the  modeling 
-  2 

error,  is  the  minimum  sum  of  error  squares  for  the  jth  model  and 

a  is  the  corresponding  parameter  vector  that  results  in  minimum  sum 


12 


of  errors  squared.  X_  is  the  design  matrix  described  in  1  . 3  •  1 « 

In  the  maximum  likelihood  procedure,  the  maximum  value  of 
the  likelihoods  from  (1.8)  are  computed  for  each  of  the  models  and 
is  taken  as  a  measure  of  the  relative  likelihood  of  the  model.  This 
approach  too  has  the  obvious  drawback  in  that  it  does  not  take  into 
account  any  asymmetry  in  the  likelihood  function.  Besides,  the 
prior  information  about  the  individual  model  form  is  generally  not 
utilized  in  the  maximum  likelihood  procedure. 

The  Bayesian  procedure  is  to  compute  the  expected  likeli¬ 
hoods  of  the  data  for  a  particular  model  form,  over  the  entire 
parameter  space  a_.  . 


L(yJ<l>j) 


L(yJ<f>j  ,  a_j  ) .  f  (a_. )  .da_j 


(1.9) 


The  symbol  L  (^J  <j> . )  is  used  to  denote  the  expected  value  of  the 
random  variable  L  (y_|  j  ,  a_. ) .  In  obtaining  the  expected  likelihood 
of  the  model,  one  may  use  the  prior  or  posterior  distribution  of  a_j  . 


There  are  arguments  for  and  against  the  use  of  the  either  distribution. 
A  discussion  on  the  subject  is  provided  in  Chapter  4. 

As  an  example,  using  the  posterior  distribution  of  the 
parameters  and  a  uniform  prior  for  the  parameters,  the  expected 
likelihood  is  [see  Plackett  (i960)] 
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Uy.Uj)  =  (Llx.ji»  aj2tIjAj 


-  (2tt)  ”n/2  .  I  x  .  C  „  .  X  .  t-4-a  . 2I I 
“J-J  “J  J 


exp .  { -  j  [rx_.a_.  ] 1  U  £  X_  t+a  2I]  ” 1  [y-^a^.  ]  }  ( 1  .  1 0) 


where 


A,  =  x/x .  and  C.  =  -^r  A.”1 
-J  -J  “J  -J  Q '  2  -J 

J 


0.11) 


X_.  being  the  n  x  m  design  matrix  for  the  jth  model,  £  is  an  n  x  n 
identity  matrix  and  |aJ  in  general  is  written  for  the  determinant 
of 

The  exponent  of  ( 1 .10)  may  be  computed  by  noting  that 
[Scheffe  (1959)] 


[y-X.a.]t[X.C.X.t+a„2l]  1  [y-X.a,] 

“  -J-J  -J-J-J  J  -  ^  “J-J 


=  1 X . C .X . t+a .2I+[y-X .a . ] [y-X . S  .]* 

'-UrJ-TJ  J..  .-J-J 


I X  .  C  .  X . t+a . 2I 
'-J-J-J  J  - 


-  1  (1.12) 


The  model  posterior  probability  can  then  be  computed  using 
Bayes  theorem 


P  j  '  •  L  (lU  j ) 

"  _  _ J _ _ 

l  pk'L(yUk) 


(1.13) 


The  expression  on  the  right  hand  side  of  (1.12)  involves 
finding  the  determinant  of  a  matrix  of  the  order  equal  to  the  number  of 
data  points.  When  the  number  of  data  points  is  greater  than  about  20, 
serious  computation  errors  arise  in  obtaining  the  determinants  for 
(1.12).  As  is  often  the  case  in  testing  equations  of  state  in  thermo¬ 
dynamics  or  certain  kinetic  models,  the  number  of  data  points  is 
usually  more  than  50  and  not  infrequently  as  large  as  200  or  more. 

To  avoid  these  computation  errors  and  large  computing  times  in 
obtaining  the  above  determinants,  an  expression  for  the  expected 
likelihood,  involving  only  the  sufficient  statistics  and  not  the 
entire  data  would  be,  indeed,  very  useful.  Such  an  expression  is 
developed  in  Chapter  2. 

1.3.3  Sequential  Design  of  Experiments 

Sometimes,  during  the  study  of  a  physical  phenomenon,  the 
data  are  not  sufficient  to  favour  any  particular  model.  In  such 
cases,  it  becomes  necessary  to  design  more  experiments  which  should 
be  carefully  planned  to  discriminate  among  the  rival  models.  The 
design  of  experiments  involves  finding  of  a  region  where  some 
particular  model  will  fit  the  data  very  well  whereas  the  other  will 
show  a  poor  fit. 

A  study  frequently  reveals  the  existence  of  such  problems 
in  which  data  in  one  region  has  negligible  information  while  in  some 
other  region,  it  might  be  substantially  more  useful  in  discriminating 
among  the  rival  models.  Box  and  Hill  (19&7)  have  given  the  example  of 
the  following  two  simple  rival  mechanisms. 
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(a)  A  +  B  -*•  C 

(b)  A->  B$C 

It  is  required  to  ascertain  one  of  the  above  two  mechanisms. 
Figure  1.2  shows  the  concentration  of  B  against  time  Xj .  Clearly, 
the  concentrations  of  B  for  large  times  are  required  to  discriminate 
between  the  two  model  forms.  Another  practical  problem  in  literature 
is  that  of  the  Water  Gas  Shift  Reaction  where  a  large  number  of 
mechanisms  have  been  considered  [Jenkins,  Kisiel,  Price  and  Rippin 
(1963)]. 

In  discriminating  among  mechanistic  models,  divergence  in 
the  responses  of  the  various  models  alone,  is  not  a  sufficient 
criterion  for  achieving  discrimination.  It  is  the  divergence  in  the 
various  responses  relative  to  the  limits  of  errors  (say  95%  confidence 
region),  that  is  important.  The  point  of  maximum  divergence  is  not 
necessarily  the  best  discriminating  point.  As  an  example,  the  light 
lines  in  Figure  1.3  represent  the  limits  of  the  errors  and  the  heavy 
lines  indicate  the  model  responses.  Obviously,  any  point  in  the 
unshaded  region  is  better  for  indicating  the  difference  in  the  models 
than  any  other  point  outside  this  region,  even  though  greater  differ¬ 
ence  in  the  responses  of  the  two  models  exist  towards  the  right  of 
the  unshaded  region. 

The  earliest  work  in  model  discrimination  was  that  of 
Williams  (1959),  Cox  (1961,  1962)  and,  Hunter  and  Reiner  (1965). 
Williams  employed  regression  analysis,  Cox  used  a  hypothesis  testing 
approach  and  Hunter  and  Reiner  took  the  designs  that  maximize  the 
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FIGURE  1.2  TWO  RIVAL  MECHANISMS  FOR  A  CHEMICAL 
REACTION  UNDER  INVESTIGATION 
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*1 


FIGURE  1.3  EFFECT  OF  ERROR  VARIANCE  IN  MODEL  DISCRIMINATION 
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deviation  between  the  predicted  values  from  the  various  models.  Sub¬ 
sequently,  Roth  (1965)  developed  a  criterion  that  calculates  the 
weighted  average  of  the  separation  from  each  given  model.  The  weights 
are  the  posterior  probabilities  of  the  models  from  all  the  prior  data 
points.  The  criterion  to  choose  the  (n+l)st  experimental  point  is  to 
maximi ze 


Z(x) 


m 

l 

i  =  1 


7T 

1=1 


Yj  (x)  -  y.  (x) 


+J 


(1.14) 


where  x_  is  the  set  of  independent  variables,  p.  n  is  the  posterior 
probability  of  the  i th  model  from  n  previous  observations  and  ^  (x) 
is  the  expected  value  of  the  prediction  from  the  jth  model  under  the 
experimental  conditions  x_. 

All  the  above  criteria  fail  to  take  into  account  uncertainty 
in  the  predicted  values  of  the  dependent  variables.  Box  and  Hill 
(1967)  developed  a  criterion  that  accounts  for  this  uncertainty  in 
the  prediction  of  the  response  of  a  model.  The  criterion  is  based 
on  the  concepts  of  information  theory,  as  developed  by  Kulback  (1959), 
and  the  idea  of  entropy  as  a  measure  of  the  randomness  or  the  uncer¬ 
tainty  of  a  situation. 

For  a  system  of  m  models,  the  entropy  at  the  end  of  n  obser¬ 
vations  may  be  defined  as 


S  = 


m 


I  p.  In  p. 

^  'in  '  1 


i  =  1 


1  ,n 


1  ,  n 


(1.15) 


The  entropy  of  a  system  is  at  maximum  if  all  the  models  are 
equally  likely  and  is  zero  if  one  of  the  models  is  certainly  true. 

Box  and  Hill  showed  that  the  negative  upper  bound  of  the 
expected  entropy  change  resulting  from  the  (n+l)st  experiment  per¬ 
formed  at  x  i s  given  by 


m 


m 


°(x)  =  7  y  i 


0  L  P  i  r>  P  5 

c  i  •  °ii  i  « n  j  « n 

i=i  j= i +i  *  j  * 


2  2 
a .  -a . 

.J ! 


,2^  2  w  2  2  v 

(c  +o .  )  (a  +a .  ) 
i  J 


(1.16) 


where  a  is  the  known  error  variance  for  a  single  measurement  and 
aj  =  ^(XX)  x^)  a  X_  i  s  defined  in  section  1.31  and  x^  is  the 

row  vector  that  would  be  the  (n+l)st  row  of  X_  for  an  observation  made 
at  x_. 

The  design  of  experiments  for  model  discrimination,  would 
now  be  to  choose  a  point  that  maximizes  D.  It  will  be  however,  more 
logical  to  have  a  criterion  that  gives  the  expected  change  in  entropy 
rather  than  an  approximation  of  it  in  terms  of  the  upper  limit  as  in 
the  case  of  Box  and  Hill  criterion.  Reilly  (1969)  obtained  an  expres¬ 
sion  for  the  expected  change  in  entropy  involving  an  infinite  integral 
that  has  to  be  solved  numerically.  The  expected  change  in  entropy 


was  given  to  be 
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m 


R(x)  =  -  7  (l  +  ln2ir)  -  £  p.  {y  ln(a2+a.2) 

•  i  I  •  I  I 
1  =  1 


+  Ej  [In  q(x,  Yn+1)> 
Yn+1 


(1.17) 


where 


yn+])  -  l  Pj ,nLi  yn+|J 

1  =  1 


L .  (x_,  y^+^  )  is  the  expected  likelihood  of  the  ith  model  for  the 

(n+l)st  observation  and  E.  [In  q(x,  y  .)]  is  the  expected  value  of 

■  —  n+ 1 

In  q  (x_,  yn  +  ]  )  over  yn+^  where  Vn+]  is  normally  distributed  with  mean 

2  2 

9n+] (x)  and  variance  c  +a.  .  Reilly  suggests  evaluating  the  last 
term  numerically  using  Gauss-Hermi te  quadrature. 

However,  it  would  be  convenient  to  have  a  simple  analytical 
expression  for  the  expected  change  in  entropy  criterion.  Such  an 
expression  is  developed  in  Chapter  3- 
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2.  EXPECTED  LIKELIHOODS  OF  DATA 


2.1 


Development  of  the  Expected  Likelihood  of  the  Data  from 

the  Sufficient  Statistics 


The  inadequacy  of  the  existing  expressions  to  compute  the 


expected  likelihood  of  the  given  data  when  a  large  number  of  obser¬ 
vations  is  involved  was  seen  in  section  1  .3=  2.  Serious  computation 
errors  arise  when  the  number  of  data  points  is  in  excess  of  about  20. 
In  addition,  the  computing  time  required  to  evaluate  the  determinants 
of  matrices  of  the  order  of  100  or  more,  can  be  considerable.  The 
difficulty  originates  from  the  fact  that  the  entire  data  set  is 
required  in  calculating  the  likelihoods  by  existing  methods.  Need 
for  a  procedure  that  would  require  only  the  sufficient  statistics  of 
the  data  is  apparent.  Such  expressions  are  developed  in  this  chapter. 
In  the  treatment  that  follows,  errors  are  assumed  to  be  normally 
distributed  with  mean  zero  and  a  known  variance.  Expressions  have 
been  developed  using  both  the  parameter  prior  and  posterior  distribu- 
t  i  ons . 

2.1.1  Expected  Likelihoods  Using  Parameter  Prior  Distribution 


When  the  errors  are  assumed  to  be  normally  distributed  with 

2 

mean  zero  and  variance  (known)  a  ,  the  joint  prior  distribution  of  the 


parameters  may  be  represented  by: 


(2.1) 
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where  m  is  the  number  of  parameters,  and  y^  are  respectively  the 
prior  variance-covariance  matrix  and  the  expected  value  of  a_,  |^| 
being  the  determinant  of  C^. 

The  likelihood  of  the  occurrence  of  n  observations  given 
the  jth  model  is  true  with  the  set  of  parameters  a_,  is 

L(y_|cj>  ,  a)  =  (2Tia2)  n/2.exp.{-  — «■} . 

J  2a 

exp  .  { -  j  [aj-a_]  t£  1  .  [a_-a_]  }  (2.2) 


where  a_  is  the  value  of  the  parameter  vector  a_,  that  gives  the  least 

A  £  XTX_ 

sum  of  errors  squared,  S  and  C_  =  — ,  X_  being  the  design  matrix  for 

o 

the  jth  model  as  defined  in  section  1.3.1. 

The  expectation  of  the  likelihood  of  the  data  over  the 
entire  parameter  space  may  now  be  written  as 


L  (yj  <f>  j  )  = 


CO 


L  (y  1 4> j  ,  a)  .  f  1  (a)  . da_ 


■  CO 


(27Ta2)"n/2.exp.{-  ^U-}.(2tt  )  m/2  .  1^1 

2a 


exp.{-  —  [a_-ja] t  .£  '  [aj-£_]  } 

exp.{-  rr  [a-y  ]  .C  ^  [a-y  ]}.da 
z  — -o  — o - o  — 


(2.3) 
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Def i ne 


Now 


S  =  [a-aj^C  ^.[a-a]  +  [a~u  ]t.C  ^  .  [a-y  1 

—  ~  — - -  —  — o  — o - o 


(2.4) 


S.  =  [a  -  [C_1  +  C’1]"1.[c‘1y  +  C_1a]]t.[c”1  +  C_1] 

—  — o  —  — o  — o  —  —  — o  — 


fa  -  [C*1  +  C'r'fc'V  +  C  ” 1  a  ]  ] 

—  *-o  —  — o  — o  —  — 


(2.5) 


and 


S0  =  atC_1a  +  y  tC_1y  -  [C  'y  +  C “ 1  a ] r  [ C " 1  +  c”1]"1 
2  - - —  —oo—o  — o  — o  —  —  — o  — 


[C  V  +  C  '  a] 

—o—o  —  — 


(2.6) 


S  =  [a-a]tC  !  [a-a]  +  [a-y  ]LC  '  [a-y  ] 

- —  — —  - o  -o - o 


t  - 


t  - 1  t  -L  „t  - 1  .  „t  - 1  „ 

=  a  C  a  -  a  C  a-a  C  a  +  a  C  a 


t  -1  t  -1  r-l  ,  t  -I 
+  aC  a-aC  y-yC  a  +  y  C  y 

—  — o  —  —  — o  — o  ~~o~~o  —  — o  o  — o 


(2.7) 


Adding  and  subtracting  [C  V  +  C  ^a]t[C  *  +  C  ^ [C  V  +  C  ^a] 

**  —^3  — o  —  —  — o  —  — n  — o  —  — 


— O 


from  the  right  hand  side 


S  = 


■c''p  +  c"'a]t[c'1 

■— o  -HD  —  —  — O 


tr-l 
a  C  a 


t„- 


+  C  “  ]  ”  1  [  c  "  1  y 

+  c"'s] 

—  — o  — o 

t  -1  t  -1 

t  - 1 

C  a  +  a  C  a 

-  a  C 
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-  y  C  a} 
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(2.8) 
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A1  so 


S  =  [a.  -  tc'1  +  c',]''[c‘ly  +  c'la]]t[c'1  +  C_l] 

[a  -  [C_1  +  C~ ' ] ~ 1 [C~ 1 y  +  C_,a]] 

—  — o  —  — o  — o  —  — 

=  [a  -  [C-1  +  c'1]"1.[c"1y  +  C  ^a]]t]]C  ^  +  C~]]a 

—  — o  —  — o  — o  — o  —  — o  —  — 


-  tcl’u  +  C~'a]] 
—o—o  —  — 
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a  [C  +  C  Ja  -  C  u  +  C  a]  a 
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+  [C  V  +  C  ^  a ] ^ [ C  *  +  C  M  * [C  V  +  C  'a] 
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t  rr"  I  ,  p - 1  ^ I 
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(2.9) 


This  is  the  same  expression  as  in  the  parenthesis  of  (2.8) 
so  that  now 


S  =  S,  +  S, 


(2.10) 


It  may  be  noted  that  only  contains  the  parameter  vector  as 


random 


variables;  S^,  is  a  constant  in  this  respect  and  may  be  taken  out  of 
the  kernel  in  the  expressions  of  (2.3).  Substituting  the  above  results 


in  (2.3)  and  taking  the  exponential  involving  out  of  the  integral, 
the  expected  value  of  the  likelihood  now  becomes 

:2 


2\~n/2  r  S _ /„  \  -m/2  i  1/2  r  2-, 

(2TTa  )  „exp.{-  — -*T  .  (2tt)  .  |C  |  .exp.{-  — }. 

2a  “°  2 


; 

—  00 


exp.{-  Jr  [a  ~  [C_^  +  C  ']  +  C  1a]]t[C^1  +  C 


2  —  —o  - 


—o—o  —  —  -o  — 


[a  -  [C  *  +  C  *]  ^ [C  V  +  C  *a]]}.da 
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(2.11) 
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But 


(2TT)‘m/2„  |c_1  +  c"]  I  i/2, 

1  — o  —  1 


exp.{-  J  [a_ 


[C_l  +  c”1  ]"'  [c"'y  +  (f'a]]1 

— o  —  — o  — o  —  — 


[C'1  +  C* 1 ] [a  -  [C_i  +  C',]',[c‘,u  +  C~ 1  a] ] } 

— o  —  —  — o  —  —0—0  —  — 


.  da_ 

=  1  (2. 12) 


Therefore  the  expected  value  of  the  likelihood  is 
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10  L  +  L 


exp .  { -  ]r  [alc  ]a  +  u  tC  V  ]} 
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1  -o  — 


exp.{+  |  [C  V  +  C  1a]t[C  1  +  C  1 ^ [ C  ^ y  +  C  ^ a] } 
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(2.13) 


The  above  expression  reduces  to  zero  for  a  uniform  improper 


prior. 
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2.1.2  Expected  Likelihoods  Using  the  Parameter  Posterior 
Di stribution 


The  posterior  distribution  of  parameters  is  given  by 


f 11  (a)  =  f  (m j  (a  |  [C~ 1  +  C~  1  ] ” 1  [C_  V  +  c” 1  a]  ,  [C~5  +  c"1]"1) 

—  N  — O  —  — O  — O  —  —  — O  — 


=  (2,a2)-m/2.|C-'  +  C-y/2 

“ O  —  1 


exp  o  { -  ]r  [a  -  [C  *  +  C  M  1  [C  \  +  C  1a]]t[C  '  +  C  '] 

7  —  — o  —  -o—o  —  _  L_o  — 


[a  -  [C"1  +  C]]_1[C"V  +  C_  1  a  ]  ]  } 

—  — o  —  — o  — o  —  — 


(2.14) 


The  expected  likelihood  of  the  data  for  the  given  model  using  the 
posterior  joint  distribution  of  parameter  is 
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(2.15) 


Let 


S  =  [a  -  [C_!  +  C”1]"1[c"1y  +  c"1a]]t[c”1  +  C_1] 

—  — o  —  —0—0  —  —  — o  — 
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—  — o  —  —o—o  —  — 
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S.  =  [a  -[C  '  +  2(f1ri[c“1u  +  2C_1a]]t[C  1  +  2C_I] 

I  —  L— O  —  -0-0  —  —  -o  — 
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Expanding  the  right  hand  side  of  (2.17)  and  regrouping 
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(2.18) 
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The  expected  likelihood  may  now  be  written  as 
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Therefore  the  expected  likelihood  for  the  given  model  is 
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'-o  —  1 

1/2 


1 


jC-i  +  C"V/2 

(2tto2)  n^2.exp.{-  — =*} .  - ; — Ypr  .exp,{-  —  a^C.  ' a_) 

2c  |C"  +  2C"  I  2 


— o 


exp.{-  y[[C  V  +  C  1a]t[C  1  +  C  ^ [C  +  C  'a] 

2  — O  — O  —  — •  — o  —  — n  — n  —  — 


-  [C~ 1 y  +  2C"1a]t[C_1  +  2C_1 ]_1 [C_1y  +  2C_]a]]}  (2. 20) 

— o  — o  —  —  — o  —  — o  — o  —  — 


. 


. 
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For  uniform  improper  prior  =  0_  and  the  above  expression 


reduces  to 


2v-n/2  r  ,  / 1 \ m/2 

(2'Tra  )  »  exp  .  { - -y) .  (y) 

2a 


(2.21) 


This  simple  expression  gives  a  theoretical  justification 

for  the  conventional  least  squares  criterion.  When  the  number  of 

-2 

parameters  m,  is  the  same  for  all  the  competing  models,  S  is  a 

quantity  that  in  itself  uniquely  defines  the  model  performance.  The 

-2 

likelihood  is  only  a  function  of  the  S  of  the  model.  In  general, 

5  m/2 

this  quantity  should  be  corrected  by  an  additional  factor  of  (— )  , 
as  to  favour  the  model  with  fewer  parameters.  This  is  evident  from 
(2.21) 


so 


3.  SEQUENTIAL  DESIGN  OF  EXPERIMENTS  FOR 
MODEL  DISCRIMINATION 
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3 • I  Development  of  Expected  Entropy  Change  Criterion  for 

Model  Discrimination 

The  objective  here  is  to  obtain  a  criterion  that  will  choose 
the  experimental  condition  x_,  such  that  the  response  of  the  system 
will  be  expected  to  result  in  maximum  discrimination  between  the 
competing  models .  We  start  with  the  concept  of  entropy  as  the  measure 
of  uncertainty  associated  with  the  knowledge  of  a  system.  Shannon 
(19^8)  introduced  this  concept  of  entropy  in  communication  problems. 
For  a  set  of  m  competing  models,  the  negative  of  expected  entropy 
change  (or  the  expected  rate  of  transmission  of  information)  for 
the  (n+l)st  observation  is 


-R  =  entropy  at  the  input  -  expected  entropy  at  the  output 


m 


l  p  .In  p  .  -  (- 1 ) 
“  r  n  ,  i  n  1 


i  = 


n+1 


m 


^Pn+l.i  Pn+1 ,  i 


f(yn+l)dVl 

(3.0 


where  f(yn+j)  is  the  unconditional  distribution  of  the  prediction 

y  ,,  and  p  .  .  is  the  resultant  probability  of  the  « th  model  after 
’ n+1  rn+l ,i 

the  (n+l)st  observation  has  been  utilized. 

After  observing  yp+^ ,  the  posterior  model  probabilities 
Pp+I  .  can  be  computed  from  Bayes1  theorem 


pn.lL(V|R|) 

’n+1  ,  i  m 

jI1Pn,jL(VllV 


(3.2) 


31 


where  L(yp+^  |cj>„),  the  expected  likelihood  of  observing  y^_Ll  given  the 


n+1 


model  form  <f>  is  true,  is  obtained  by  integrating  the  likelihood 

L  (yp+j  I  <f>  •  »a_. )  of  observing  y  ^  i  f  <j>  „  is  true  with  parameters  a_t 
J  J  ^  J  J 

over  the  entire  parameter  space  }  that  is 


L(yn+llV  = 


L  (y  ,  j<j>  „  ,a  .  )f  (a  .  )da  . 
n+1 1  j  — j  — j  — j 


(3.3) 


and  f  (a_j )  being  the  prior  parameter  distribution  which  is  multivariate 

normal  and  incorporates  the  n  prior  data  points.  From  the  assumption 

2 

that  the  errors  are  normally  distributed  with  known  error  variance  c 


L(yn+)kj.,aj)  -  Vyn+i lvn+i ,j>°  ) 


(3.4) 


A 

where  y^+.  j  is  the  expected  value  of  the  prediction  from  the  jth 
model  at  the  value  of  £n+|  •  corresponding  to  the  experimental 

condition  where  yn+^  has  been  observed.  Substituting  (3.4)  into  (3. 3) 
leads  to 


L(yn+ll*j)  '  fN(yn+|!yn+1J.°2  + 


(3.5) 


where 


2  /vt  v  v-1  t  2  v 

o  ,  .  .  =  x  .  (X  .X  j  x  .  .c 

n+l,j  -n+1  ,j  — n  ,  J— n  ,/  -n+l,j 


(3.6) 


as  defined  in  section  1.3.3  the  term  a  accounts  for  the  error  due  to 


z 

a  single  measurement  and  a  .  describes  the  error  due  to  the 

n+!,j 


uncer¬ 


tainty  in  the  parameters  given  that  the  model  form  <{> .  is  correct;  . 
is  the  design  matrix  of  the  n  prior  observations  and  x^+^  j  is  the 
design  vector  for  the  jth  model  under  the  experimental  conditions  for 
observing  yp+^ ,  as  defined  in  sections  1.3.1  and  1.3*3. 
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After  Vn+i  has  been  observed,  all  the  quantities  L  (y  n_j_  ^ 
and  Pn+|  j  are  scalars.  However,  before  the  (n+l)st  observation  has 
actually  been  made,  we  treat  y  ^  as  a  random  variable  and  L  (y^+ j  |  .  ) 
and  PR+^  j  in  (3.2)  are  now  functions  of  a  random  variable  and  are 
themselves  random  variables.  Furthermore,  it  is  assumed  that  the  set 
of  models  under  consideration,  which  we  call  a  metamodel,  is  mutually 
exhaustive,  i ,e.  the  data  is  generated  by  a  model  within  the  metamodel. 
Of  course,  in  many  situations,  it  is  always  possible  that  some  other 
model  may  describe  the  data  better  than  all  the  models  in  metamodel. 
However,  we  cannot  take  into  account  that  which  we  do  not  perceive. 
Given  this,  it  follows  that 

m 

f(yn-H)  =  J,Pn,jL(yn+lhj)  (3‘7) 

J=1 

that  is,  the  unconditional  density  of  Yn+j  is  a  weighted  average  of 
the  likelihoods  of  generating  the  data  from  the  various  models  under 
consideration.  This  assumption  enables  the  exact  solution  of  the 
i ntegra 1  in  (3.1) 

Substituting  (3.7)  and  (3.2)  in  (3.1) 

m  r 

-  R  =  -  I  p  .In  p  „  + 

J  =  1  n ,  !  n’!  J 

Vi 

•  dyn+i  (3.8) 


m 


£pn,iL(yn+l 


Uj  ) 


pn,iL(yn+llh) 


The  integral  in  (3.8)  may  be  expanded  to  give 
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n+1 


m  m  L  (y  ,  ]  cj) .  ) 

Ip  ;  *-  (y  a.  i  I  <f>  -  )1  n  P  •  +  /  p  •  L  (y  ,  ,  U » ) »  1  n  v;-n - r— 

",  n,i  'n+lIYr  n ,  i  >,Kn,i  7n+i1Yi  f  (y  .) 


dy 


n+1 


(3.9) 


S  i  nee 


L  (yn+|  Uj )dVn+i  =  ]>  (3.9)  becomes 


y 


n+1 


m 

I  p  .In  p  .  + 

j  =  I  n  ,  i  Kn , i  , 


n+1 


m 


L(yn+I l+s> 


I  p  .  L  (y  ,  j  4> .  )  1  n  -  v-/-n- - “r 

+  1  n,i  7n+l|Yr  f  Cyn+1) 


dyn+,  (3.10) 


Substituting  in  (3.8) 


-  R  = 


m 

I  p 

i=i 


n,  i 


L(V||*1) 


L(Vll*l)1n"f(yn+1)  dy 


y 


n+1 


n+1 


(3.H) 


From  (3.6)  and  (3.7),  the  unconditional  distribution  of  the 

-  2 

dependent  variable  yn+^  may  be  written  in  terms  of  yn+j  and  an+]  where 


m 


yn+l  .  L ,  Pn  ,iyn+l  , 
i  =  l 


(3.12) 


and 


2  _  r  2 

°n+l  •  ,i°n+l  , 


(3.13) 


So  that 


f(Vl>  =  VVl'Vl’  °2  +  °n+l } 


(3.14) 


3*t 


Utilizing  the  above  relations,  the  integral  in  (3.11) 

reduces  to 


,  t  lx  m  L(Vn+l 
L  yn+l  I  ^  j }  n  T^") . -  d>k+l 

f 

=  ,  L(yn+lK)  ,n  L(yn+dh)dVl 

yn+l 

/ 

",  L(Vl  I*;)  ln  f(yn+,l*0)dVl  (3J5) 

yn+] 


The  first  integral  on  the  right  hand  side  of  (3.15)  can  be  evaluated 
as  fol lows 


L(*n+ii*i)  { 


i  I n2u (a2  +  o^+]>.) 


n+1 


1  yn+l  ■  yn+l , I >  ,, 

2  ,'2,2 - ~~ }dyn+l 

(0  +Vi,i) 


/  A  \  Z. 

The  expected  value  of  (y^  “  Yn+^  j)  with  respect  to  Yn+p  assuming 

,2  2 

model  i  to  be  correct,  i s  (a  +  °n+^  j)°  Consequently,  the  above 
integral  becomes 


k  [l  +  ln2l,(°2  +  °n+l,i)] 


(3.16) 


The  second  integral  on  the  right  of  (3.15)  can  be  written 


as 


L(yn+1  k;  )  ”  ln27T  (c2  +  a2^, ) 


1  (yn+l  "  yn+l ^ 


n+1  2 


n+1 


■■y  ;■ "  2— 

°  °n+ 1 


-}dy 


n+1 


(3.17) 


■ 
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A'so,  (yn+,  -  vn+|)  -  (yn+)  -  vn+lji)  +  (yn+l>i  -  vn+l) 


+  2 (yn+ 1  ‘  Vl,l)(Vl  '  yn+l) 


(3.18) 


The  expected  value  of  the  last  term  in  (3.16)  is  zero  and 
the  second  term  is  just  a  constant.  (3.17)  now  simplifies  to 


.  .  (_2  2  ,  0  °n+l , i  (yn+l,i  "  Vl} 

lnMo  +  °n+l  +  2  "  2 - +  - 2 - 2 - 


a  +  a 


n+1 


a  +  a 


n+1 


(3.19) 


Substituting  the  results  of  (3.16)  and  (3-19)  Sn  (3.8) 


m 

R  =  0.5  Ip 

u  r  n 

i  =  l  ’ 


2  2  2  2 

a  +  a  ,  ,  o  a  ,  ,  .  -  a  ,  , 

,  n-H  1,1  n+1  ,  i  n+1 

-  1  n  — -  ■  „  ■-  +  — ^ — x — -  + 


(y 


n+1  , 


yn+l> 


2  2 
a  +  an+l 


2  2 
0  +  °n+l 


2  ^  2 
a  +  Vi 


(3.20) 

Now,  the  second  term  of  (3.20)  disappears  as  can  be  seen  from  the 
relation  (3.13).  Thus 


-  R 


m 

0.5  l  P 
i  =  1 


n ,  i 


2  2 
a  +  a  ,  .  . 

n+1 , i  . 

n  — 5 - ”2  + 

°  +  °n+ 1 


(y 


n+1  ,  i 


yn+l 


2  2 
a  +  an+1 


(3.21) 


The  expression  appeals  in  its  simplicity.  It  also  supports 

the  attempt  of  earlier  workers  to  maximize  the  difference  in  the 

2 

response  of  various  models  by  the  inclusion  of  (y  .  .  -  y  ,)  in 

n+ I , i  n+1 

2  2 

(3-21 )  but  restrains  it  with  the  uncertainty  term  (a  +  a  ,)  It  is 

n+ 1 

the  ratio  of  these  two  terms  that  is  important.  In  most  cases  studied, 
the  first  term  of  (3.21)  was  significantly  less  important. 
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y 


n+1 


When  one  model  is  known  for  certain,  c  .=  a  ,,  .  and 

*  n+1  n+1 , i 

A 

y  „,  where  the  ith  model  is  the  correct  one.  Also  p  .  =  0 
n+ l,i  n ,  j 

this  gives 


R  =  0 


which  is  expected. 

For  comparison,  a  restrospect i ve  view  of  the  two  most 
important  earlier  design  criteria  due  to  Roth  and  Box  and  Hill  may 
be  provided  here.  Consistent  to  the  notation  used  here,  Roth's 
criterion  is  to  choose  the  (n+1 ) s t  experimental  setting  that  maximizes 


m 


z(x)  =  l 


i  =  1 


n,  i 


m 

IT 


j  +  1 


|y 


n+1 , j  - 


(x)  - 


y 


n+1  ,  i 


(x)  | 


Box  and  Hill  required  maximizing  the  negative  upper  bound 
of  the  expected  entropy  change 


D(x)=l  I  l 

i  =  l  j  =  i  + 1 


m  in 


p  .  p 

n, i  n  ,  j 


i 


2  _  2 
°n+l , i  °n+i , j 


(a2  +  c2  ,  . )  (a2  +  a2  .) 
n+1 , i  n+1  ,  j 


Vi.i'i* 


yn+l , j 


1 


2  2 
o  +  a  .  . 

n+1  ,  i 


2^2 
a  +  a  , ,  . 

n+1 ,xl 


It  may  be  noted  that  Roth's  criterion  failed  to  take  into 
account  the  uncertainty  in  the  response  and  was  based  on  more  or  less 
an  empirical  interpretation  of  model  discrimination.  Box  and  Hill's 
criterion  was  developed  from  the  concept  of  expected  entropy 
change  but  is  not  exact.  An  upper  bound  of  the  theoretically 


desirable  ‘expected  entropy  change'  is  used  as  an  index  for  the 
design  of  experiments,.  The  proposed  criterion  in  (1.21)  is  both 
simple  and  exact  without  involving  any  further  restrictions. 
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3 • 2  Examples 

The  examples  chosen  to  illustrate  the  use  of  expected 
change  in  entropy  criterion,  have  all  appeared  in  literature 
[Hill  (1966)].  The  first  two  examples  involve  linear  models  and 
the  last  one  consists  of  discriminating  between  non- 1 i nea r  kinetic 
model s . 


Examp  1 e  3 ■ 1 


Consider  discriminating  between  the  following  polynomial 


mod  els. 


Model 

1  . 

v0) 

.  a(l) 
al 

Model 

2. 

y(2) 

.  a  (2) 
d| 

Model 

3» 

y(3) 

.  a  (3) 
a] 

Mode  1 

4. 

yW 

.  ,(« 
al 

x 


♦  a<2)x 

+  a‘3)x  ♦  a^x2 
x  +  a^x2 


The  above  four  models  represent  respectively,  a  straight 
line  through  the  origin,  a  straight  line  not  necessarily  through  the 
origin,  a  general  quadratic  and  a  quadratic  forced  through  the  origin. 


The  data  was  generated  from  model  3  with  a,w/  =  =  a^~"  =  1  and 


(3)  _  (3)  =  (3) 

2  3 


error  variance  of  a  =1.  The  independent  variable  was  constrained 

to  0  <  x  <  10  and  integer  values.  Starting  with  equal  probabilities 

for  all  the  models,  that  isp  ,  =  p  „  =  p  ->  =  P  \  ~  '25  and 

o,l  o,2  o,3  o,4 
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using  five  data  points  at  x  -  0,  1,  2,  3  and  4,  the  model  probabilities 

-  c  _c 

were  p  .  =  2.5  x  10  ,  p  9  =  6.9  x  10  ,  p  -  =  .445  and  p  ,  =  .554. 
The  experimental  space  was  then  scanned  to  find  out  the  value  of  the 
independent  variable  that  gives  maximum  expected  entropy  change. 

The  results  are  given  in  Table  3-1.  Actual  change  in  entropy  was 
calculated  to  compare  with  the  expected  change  in  entropy.  The 
results  given  here  are  typical  of  various  runs.  Total  actual 
entropy  change  was  generally  found  close  to  the  total  expected 
change  in  entropy  for  the  entire  run. 

In  this  particular  example,  the  optimum  experimental  point 
was  always  found  to  be  zero,  which  is  clearly  expected  from  the  model 
forms  of  the  two  important  competing  models,  that  is  model  3  and 
model  4.  The  posterior  density  of  a_ was  used  in  computing  the  model 
probab i 1 i t i es . 

Example  3»2 

This  example  is  to  illustrate  the  mode!  discrimination 
procedure  when  some  more  complex  models  are  reducible  to  the  true 


model  form. 
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Table  3 • 1 

Results  When  Discriminating  Among 
the  Four  Polynomial  Models 


(Mode  1 


n 

X 

pn,l 

Pn  9 
n  ,2 

0 

.25 

.25 

5 

- 

.001 

6 

0 

- 

.001 

7 

0 

- 

.001 

8 

0 

- 

.001 

9 

0 

- 

.001 

10 

0 

- 

.001 

1 1 

0 

- 

.001 

12 

0 

_ 

.001 

3  True) 


Pn,3 

Expected 
Entropy  Change 

P  i  "R 

Kn  ,4 

Actua 1 
Entropy 
Change 

.25 

.25 

.445 

.554 

.881 

.118 

.151 

.272 

.963 

.036 

.272 

.209 

.948 

.051 

.202 

1 

O 

OO 

.943 

.056 

.037 

-.013 

.942 

.057 

.029 

-.003 

.962 

.037 

.024 

.062 

.997 

.002 

.016 

.142 

bO 


Taking  the  same  set  of  models  as  in  Example  3*1,  it  is  clear  that 

(3) 

model  3  is  reducible  to  model  2  by  putting  a^  =  0.  Thus  if  the 

data  were  generated  by  model  2,  model  3  will  fit  the  data  equally 

well  and  it  might  be  thought  impossible  to  discriminate  between 

these  two  model  forms.  In  fact,  the  procedure  does  discriminate 

between  the  two  models  to  favour  the  model  with  the  fewer  number  of 

parameters,  the  reason  being  that  the  average  variance  of  the 

prediction  from  the  model  with  the  least  number  of  parameters  is 

smallest.  It  is  also  evident  from  our  new  expressions  <^f  the 

1  2~ 

expected  likelihood  in  Chapter  2,  where  a  factor  of  (— )  gives  the 

maximum  expected  likelihood  for  the  model  with  least  m,  when  all 

-2 

such  models  have  the  same  S  . 

The  data,  in  this  example,  was  generated  from  model  2 

(2)  (2) 

with  a^  =  a^=  1  and  unit  error  variance.  To  start  with,  all 

models  were  considered  equally  likely  and  first  5  points  data  gave 

the  posterior  model  probabilities  of  p  =  .32,  p  ?  =  .251, 

5>  1  5>z 

pc  _  =  .187  and  p_  ,  =  .233.  The  results  are  shown  in  Table  3-2. 
b  j  3  5  •  ^ 

There  is  a  marked  preference  for  the  simpler  model  (model  2)  after 
12  additional  data  points,  although  more  runs  are  required  for 


further  discrimination. 


Table  3*2 
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Results  When  Discriminating  Among  the  Four  Polynomial  Models 
(Model  2  True,  Model  3  is  Reducible  to  Model  2) 


Expected  Actual 

n  x  Pn  1  pn  2  Pn  3  Pn  4  Entropy  Change  Entropy  Change 


0 

.25 

.25 

.25 

.25 

- 

5 

.329 

.251 

.187 

.233 

- 

6 

0 

.179 

.300 

.163 

.358 

.70 

.033 

7 

0 

.304 

.377 

.109 

.209 

.13 

.033 

8 

0 

.038 

.742 

.194 

.026 

.113 

.540 

9 

0 

.003 

.775 

.221 

.002 

.052 

.20 

10 

0 

.002 

.759 

.237 

.002 

.005 

-.016 

1  1 

0 

.003 

.737 

.258 

.002 

.004 

-.029 

12 

0 

- 

.733 

.267 

- 

.005 

.022 

13 

4 

- 

.745 

.255 

- 

.001 

.012 

14 

10 

- 

.779 

.221 

- 

.001 

.039 

15 

4 

- 

.778 

.222 

- 

.004 

-,001 

16 

4 

- 

.750 

.249 

- 

.004 

-.032 

17 

4 

- 

.805 

.195 

- 

.006 

.068 

' 

Example  3-3 


Suppose  an  investigator  is  studying  a  set  of  non-linear 
models,  it  could  be,  say,  the  four  different  mechanisms  of  the 
reaction  A  ->  B,  depending  upon  whether  the  reaction  is  first, 
second,  third  or  fourth  order.  The  corresponding  models  to 
describe  these  mechanisms  shall  be 


Mode  1 

1  . 

y(,)  - 

exp .  { 

-  a  j1  ^  exp  (-  a^1  Vx2)  } 

Mode  1 

2. 

II 

CM 

{1  + 

a  |2^exp.  (-  a22  Vx2) }  1 

Mode  1 

3. 

-< 

OJ 

II 

{1  + 

2a^exp.(-  a^/x2)r1/2 

Model 

4. 

(4) 

y  - 

{1  + 

3a1(3)exP.(-  a2(%)f1/3 

s  the 

react 

ion  time, 

x  is 

the  temperature,  y^'^  is 

the 

concentration  of  the  reactant  and  aj  and  a^  are  the  parameters 
of  the  models.  It  is  desired  to  determine  the  most  likely  model  by 
the  Bayesian  model  discrimination  procedure.  The  independent 
variables  are  constrained  by  the  process  to 


0  <_  Xj  <_  1 50  and 
450  <_  x2  <_  600 

A  grid  spacing  of  25  is  suggested  for  scanning  the  experi¬ 
mental  space.  The  parameter  values  are  known  to  be  between 

300  <_  a  | ^  <_  400  and  5000  <_  a^'  ^  6000  for  all  i  =  1,  2,  3  and  4. 

2 

Variance  of  the  experimental  measurements  a  is  .0025.  In  this 
example,  two  cases  were  investigated;  (a)  model  2  is  correct,  (b) 


model  3  is  correct. 
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Both  times,  data  was  generated  with  a 


f 1 ^  =  350  and  a^' }  =  5500. 


All  models  were  taken  to  be  equally  likely  at  the  start  and  a  prior 
parameter  variance-covariance  matrix  was  taken  to  be 


4 

10  0 
0  106 

A  scheme  due  to  Marquardt  (1963)  was  utilized  to  obtain 

A 

a_,  the  least  square  estimators  of  the  parameters. 

Table  3.3  and  3.4  show  the  results  of  cases  (a)  and  (b) 
respectively.  The  discrimination  was  very  rapid  and  favoured 
decisively,  the  right  models  within  3  observations.  In  both  the 
cases,  the  exploration  for  the  first  experimental  point  gives 
maximum  expected  entropy  change  at  (100,  600)  and  thereafter  the 
cornermost  point  (150,  600)  is  favoured.  To  illustrate  the  extent 
of  discrimination  in  this  example,  Table  3*5  gives  the  expected 
entropy  change  for  the  various  grid  points  for  the  search  of  second 
point  in  case  (b).  The  results  show  sharp  discrimination  at  the 
higher  values  of  x^  and  x^  and  no  d i scr imi nat i on  at  all  when 
x^  =  0.  In  another  example  (model  2  true)  40  data  points  were  gen¬ 
erated  at  (150,450)  and  the  final  results  were  found  to  be  p^g  ^  =  .56, 

p40  2  =  ,Zt1,  p4o  3  =  ,0^  and  p40  4  =  °*  This  illustrates  the 
tremendous  difference  in  the  information  contained  in  the  experiments 

that  are  carried  out  at  (150,  450)  and  those  at  (150,  600) .  The 
utility  of  the  correct  experimental  design  is  obvious  in  such  situations. 


By  comparing  the  results  obtained  by  Hill  (1966)  for  this 
example,  it  was  found  that  the  design  point  given  by  the  Box  and 
Hill  criterion  and  that  given  by  the  expected  entropy  change  criterion 
proposed  here,  are  not  always  the  same. 
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Table  3-3 

Results  When  Discriminating  Between  the 
Four  Kinetic  Models  (Model  2  True) 


n 

xi 

x2 

Pn  i 
n  9  1 

Pn,2 

Pn,3 

pn,4 

0 

.25 

.25 

.25 

.25 

1 

100 

600 

.001 

.588 

.390 

,021 

2 

150 

600 

- 

.588 

.390 

- 

3 

150 

600 

- 

.999 

,001 

- 

k 

150 

600 

- 

1  .000 

- 

- 
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Table  3-4 


Resu 1 ts 

When 

D i  scr imi nat i ng 

Between  the 

Four 

Kinetic  Models  (Model 

3  True) 

X1 

X2 

Pn  1 
n ,  1 

Pn  9 
n  ,2 

Pn,3 

Pn,i, 

.25 

.25 

.25 

.25 

100 

600 

io"6 

,044 

.449 

vn 

o 

oo 

150 

600 

- 

CO 

r-''. 

o 

.855 

.067 

150 

600 

- 

,007 

.989 

.003 

3 


X1 

0 

50 

150 

0 

150 

0 

150 

0 

50 

100 

150 
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Table  3-5 

Search  for  the  Best  Experimental  Conditions 
While  Discriminating  Between  the  Four  Kinetic 
Models  (Model  3  True)  Expected  Change  in 
Entropy  at  Various  Operable  Conditions 


x 


2 


Expected  Change 
in  Entropy 


450 

.0 

450 

CO 

o 

o 

450 

.018 

500 

.0 

500 

.021 

550 

.0 

550 

.207 

600 

.0 

600 

.156 

600 

.350 

600 

.451 
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4.  STANDARDIZATION  OF  THE 
BAYESIAN  MODEL  DISCRIMINATION  PROCEDURE 

4 . 1  Examination  of  the  Bayesian  Procedure 

The  Bayesian  method  of  discriminating  among  the  rival  models 
by  using  the  expected  likelihood  of  the  data  to  compute  the  posterior 
model  probabilities  has  been  recognized  as  the  best  available  technique 
for  model  discrimination.  By  using  the  expected  likelihood  rather 
than  the  maximum  likelihood  as  the  criterion  for  discrimination,  any 
asymmetry  in  the  likelihood  function  is  accounted  for.  Furthermore, 
unstructured  information  in  the  form  of  prior  beliefs  about  the  model 
form  or  the  parameters  can  be  utilized  by  assigning  a  prior  distribution 
to  the  parameters  of  each  competing  model  [see  Quon  (1970)]. 

However,  some  areas  of  the  Bayesian  procedure  are  not  well 
defined  and  are  controvers i a  1 .  This  chapter  is  written  with  a  view 
of  rationalizing  the  Bayesian  model  discrimination  procedure  to  some 
extent . 

A  few  points  worth  analysing  are:  (a)  assumptions  associated 
with  the  errors  and  estimates  of  error  variance  used  in  computing  the 
likelihoods,  (b)  whether  prior  distribution  or  the  posterior  distri¬ 
bution  of  parameters  should  be  used  in  computing  the  expected  likeli¬ 
hoods,  (c)  will  a  given  amount  of  information  containing  various  sets 
of  data,  always  lead  to  the  same  level  of  discrimination,  regardless 
of  the  order  or  the  groupings  of  data,  when  used  piecemeal  to  discrimi¬ 
nate  the  models  sequentially.  These  problems  are  investigated  here 
and  a  rationalized  procedure  is  recommended  in  the  next  section. 


4.1.1 


Errors  and  Error  Variance 


Suppose  the  model  under  consideration  is 

y  =  <t>  U,  a)  +  e 

where  e  represents  the  error  between  the  actual  observed  values  of 
the  dependent  variable  and  the  one  predicted  by  the  model.  This 
error  term,  then  includes  experimental  errors  and  the  errors  due 
to  the  predictive  model  being  not  the  true  model.  The  first  kind 
of  errors,  the  experimental  errors,  can  never  be  entirely  eliminated 
whatever  the  precision  of  the  experimental  measurements.  The  modeling 
error  shall  also  be  almost  invariably  present,  since  the  mathematical 
models  however  sophisticated,  can  rarely  describe  a  complex  physical 
phenomenon  exactly.  The  importance  of  understanding  the  nature  of 
these  errors  is  obvious,  since  model  discrimination  itself  is  the 
statistical  processing  of  the  errors.  No  discrimination  problem 
would  be  present  if  errors  were  not  involved. 

The  errors  are  conveniently  assumed  to  be  normally  distri¬ 
buted  in  most  of  the  discrimination  problems.  This  is  not  a  very  far¬ 
fetched  assumption  as  far  as  the  experimental  errors  are  concerned. 

"Everybody  believes  in  the  law  of  errors,  experi¬ 
menters  because  they  think  it  is  a  mathematical 
theorem,  the  mathematicians,  because  they  think  it 
is  an  experimental  fact." 

However,  there  does  not  seem  to  be  any  justification  for 
treating  the  modeling  errors  to  be  normally  distributed.  In  fact, 
it  might  be  more  apt  to  describe  the  errors  due  to  model  form  by 
simpler  functions  like  polynomials.  This  is  the  basis  of  model  improvement 


. 

I 
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in  generating  spline  functions. 

Reconsidering  the  question  of  modeling  errors,  it  might  be 
argued  that  if  modeling  errors  were  small  as  compared  to  the  experi¬ 
mental  errors,  there  probably  would  not  be  much  harm  in  assigning 
a  normal  distribution  to  the  overall  errors.  The  modeling  errors 
would  be  negligible  in  case  of  a  model  describing  the  mechanism 
reasonably  well  and  would  be  relatively  small  for  all  the  competitive 
models.  Any  model  form  that  has  large  inherent  modeling  errors 
associated  with  it,  would  be  reduced  to  a  non-competitive  position 
outright  and  thus  would  not  affect  further  discrimination  between 
the  crucial  models. 

Besides  the  advantages  of  treating  the  errors  to  be  normal¬ 
ly  distributed  are  great  because  of  the  analytical  tractability  of 
the  normal  form.  Almost  the  entire  present  discrimination  literature 
has  been  developed  starting  with  this  assumption.  Besides,  if  the 
normality  assumption  is  not  used,  then  there  arises  the  difficulty  of 
specifying  an  alternative  distribution.  The  computational  disadvan¬ 
tages  far  outweigh  any  advantages  that  might  accrue  in  terms  of  better 
discrimination.  It  is  clear  from  the  foregoing  discussion  that  it  is 
best  to  live  with  the  assumption  of  the  normality  of  errors. 

However,  it  might  be  worthwhile  examining  the  actual  errors 
systematically  over  the  entire  operating  region  for  any  correlation 
between  the  errors  and  the  independent  variables.  In  general,  correl¬ 
ated  errors  signify  the  existence  of  large  modeling  errors.  Of  course, 
the  errors  could  also  be  sometimes  correlated  due  to  a  poor  experimen¬ 
tal  technique,  but  these  can  often  be  explained  by  examining  the 
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various  measurements  and  experimental  conditions.  In  such  cases, 
when  the  model  is  found  to  be  the  cause  of  correlated  errors,  the 
model  can  be  further  improved  by  the  addition  and  or  modification 
of  the  functions.  This  can  be  particularly  valuable,  where  the 
process  is  poorly  understood  or  very  precise  data  is  available. 

It  is  to  be  noted  that  the  modeling  errors  considered 
above  are  only  errors  due  to  the  model  form  and  do  not  include  the 
errors  due  to  uncertainty  in  the  parameters  of  a  particular  model 
form.  The  errors  in  the  prediction  after  (n-l)  observations,  due  to 
the  uncertain  parameters  are  also  taken  to  be  normally  distributed 
with  mean  zero  and  variance 


2 

c 

n,  i 


t/vtvx-l  2 
x  (XX)  x  a 

— n - — n 


as  described  in  section  1.3.3-  The  normality  of  the  errors  due  to 
the  uncertainty  about  the  parameters,  directly  follows  from  the 
previous  assumption  of  normality  about  the  other  two  types  of  errors. 

Having  decided  to  treat  all  errors  as  normally  distributed, 
we  now  examine  the  magnitude  of  these  errors.  The  mean  of  these 
errors  is  taken  to  be  zero.  The  variance  of  the  errors  may  be  known 
or  unknown,  although,  in  general  it  is  very  hard  to  have  a  reasonable 
estimate  of  the  error  variance.  While,  the  experimental  errors  are 
amenable  to  some  kind  of  assessment,  the  extent  of  modeling  errors  is 
almost  never  known. 

There  are  standard  methods  of  analysis  both  when  error  vari¬ 
ance  is  known  or  unknown  [Raiffa  and  Schlaifer  (1961),  Lindley  (1965), 
Quon  (1970)].  However,  the  popular  method  is  to  obtain  an  estimate 
of  the  error  variance  from  the  least  sum  of  errors  squared  and  treat 
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it  is  known  in  computing  the  parameter  posterior  distribution  and  the 
expected  likelihoods. 

Three  types  of  estimates  of  error  variance  are  commonly 
employed,  viz.  (a)  A  prior  estimate  obtained  from  estimating  the 
measurement  errors  and  neglecting  the  contribution  due  to  the  modeling 
errors.  (b)  Maximum  likelihood  estimator  of  the  error  variance  for  a 
given  model,  computed  as  the  least  sum  of  errors  squared  divided  by 
number  of  data  points,  that  is  S  /n.  This  can  be  interpreted  as  the 
average  least  errors  squared  per  observation  and  is  itself  an  index 
of  the  model  performance.  (c)  The  unbiased  estimator,  computed  as  the 
least  sum  of  errors  squared  divided  by  the  number  of  observations  less 
the  degrees  of  freedom  associated  with  the  model,  i.e.  S  / (n-m)  where 
m  is  the  number  of  parameters  of  the  model. 

The  relative  utility  of  these  error  variances  was  examined 
for  a  few  polynomial  model  systems.  It  was  found  that  using  the 
updated  unbiased  estimate  of  the  error  variance  resulted  in  most 
rapid  discrimination  followed  closely  by  the  results  obtained  from 
the  maximum  likelihoods  estimate  of  the  error  variance.  Using  the 
error  variance  that  was  used  to  generate  the  data  gave  a  relatively 
poor  discrimination.  An  illustrative  example  is  given  here. 

Examp  1 e  4.1 

The  following  two  simple  models  were  used  for  this  example. 


Model  1 . 


(1) 


(2) 


+ 


+ 


x 

x  +  a 


(2)  2 

3 


Model  2. 
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Table  4.1 

Results  When  Discriminating  Between  the  Two  Polynomial  Models 
Using  Various  Estimates  of  Error  Variance 
(Model  2  True) 


n 

? 

Probability  of  Model  2  When  c  = 

1 

S  /n 

S2/(n-m) 

0 

.5 

.5 

.5 

5  (Prior) 

.542 

.542 

.542 

6 

.589 

.632 

.644 

7 

.642 

-^1 

00 

V£> 

.811 

8 

.687 

.892 

.913 

9 

.  808 

.948 

.963 

10 

.911 

vO 

cr> 

.994 
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(2) 

Data  was  generated  from  model  2  with  =  a„  =  a„  =  1 


(2)  _  (2) 
'2  3 


and  error  variance  of  unity.  The  posterior  distribution  of  the 
parameters  was  used  to  compute  the  likelihoods.  Updated  error  variance 
has  been  used  here.  Table  4.1  shows  the  results. 


4.1.2  Parameter  Joint  Distribution 

What  parameter  distribution,  prior  or  the  posterior,  should 
be  used  in  computing  the  expected  likelihoods?  The  question  has  been 
debated  frequently. 

Since  the  expression  for  the  likelihood  contains  all  of  the 
posterior  information,  it  seems  justified  to  use  only  the  prior  infor¬ 
mation  in  the  parameters  to  compute  the  expected  likelihood  by 
integrating  the  likelihood  expression  with  respect  to  the  prior  dis¬ 
tribution  of  the  parameters.  This  might  even  sound  analogous  to  the 
statement  of  Bayes1  theorem  that  the  posterior  probability  of  a 
hypothesis  is  proportional  to  the  product  of  prior  belief  about  the 
hypothesis  and  the  posterior  likelihood  of  the  data  given  the  hypothesis 
is  true. 

On  closely  examing  the  problem,  it  is  at  once  obvious  that 
the  Bayes  theorem  cannot  be  applied  here.  The  posterior  expected 
likelihood  of  the  data  is  simply  the  average  of  the  posterior  likelihood 
function  over  the  entire  parameter  space.  The  correctness  of  the 
results  will  depend  upon  how  good  is  the  estimate  of  the  parameter 
distribution  employed  in  calculating  such  an  average.  Clearly,  the 
posterior  distribution  of  the  parameters  is  the  best  estimate  we 
have  to  describe  the  parameter  space.  Thus,  using  the  posterior 
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distribution  of  the  parameter  sounds  more  logical  and  is  not  equiv¬ 
alent  to  'using  the  information  twice'. 

Further,  the  use  of  prior  parameter  distribution  fails  to 
give  any  results  for  an  improper  uniform  prior  since  the  expected 
likelihoods  obtained  are  zero  for  all  the  models.  Also,  it  was 
suspected  that  for  the  vague  priors,  the  model  discrimination  would 
be  very  sensitive  to  the  degree  of  vagueness  if  the  parameter  prior 
distribution  is  used.  The  following  two  examples  were  constructed 
to  demonstrate  the  effect  of  vague  priors. 

Example  A. 2 

Model  1 .  y  =  a  ^  +  b  x 

Model  2.  y  =  +  b 

1  2 

Data  were  generated  from  y  =  10  +  10  x  +  e  with  known  error  variance 
2 

c^  =  1;  5  points  were  obtained  at  x  =  0,  0.5,  1.0,  1.5  and  2.0. 

Prior  parameter  distributions  for  the  rival  models  were 

4 

10  0 

4 

0  10 

2  was  varied  to  examine  the  effect  on  the  posterior  probabilities 
of  the  models.  The  prior  probability  of  each  model  was  taken  to  be  0.5. 

The  posterior  probability  of  model  1  varied  from  .42  to  .96 
as  y  ?  was  given  values  from  [10,  10]  to  [200,  200]  when  using  the 
parameter  prior  distribution.  On  the  other  hand,  using  the  parameter 
posterior  distribution,  there  was  practically  no  effect  of  different 


u  =  [10,  10]  ; 
-o,  1 
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Table  A  .2 

On  the  Model 

Effect  of  Prior  Parameter  Values 

Discrimination  Between  Two  Rival  Polynomial  Models 

^o,2 

Posterior  Probabilities  of  Model  1  Using 

Parameter  Parameter 

Prior  Distribution  Posterior  Distribution 

[  10  10] 

.420  .261 

[  20  20] 

.428  .261 

[  50  50] 

.467  .261 

[100  100] 

.630  .261 

[200  200] 

.965  .262 

.262 
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Table  k.3 

Effect  of  Prior  Variance-Covariance  Matrix 
On  the  Model  Discrimination 
Between  Two  Rival  Polynomial  Models 


Posterior  Probabilities  of  Model  1  Using 

Parameter  Parameter 

Prior  Distribution  Posterior  Distribution 


.106 


.365 


103  0 
0  103 
102  0 
0  102 


.5M 


.  924 


.365 


.366 
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degrees  of  vagueness  on  the  posterior  probabilities.  The  results  are 
presented  in  Table  4.2. 


Example  4.3 

The  same  competing  models  were  used  as  in  Example  4.2  but 

2 

data  were  generated  from  model  1  with  a ^  =  b ^  =  1  and  a ^  =  1  ,  at 
x  =  0,  0.5,  1.0,  1.5  and  2.0.  The  model  parameter  priors  are 


U  i  =  ^  9  =  H  ,  1  ]  ancl  C  9 
— O  ,  I  ~0 , 2  — o , 2 


10- 


0 


0 


10- 


c 

— o , 


io2 

0 

io4 

0 

was  varied  from 

to 

0 

io2 

0 

IO4 

and 


was  found  to  reduce  the  model  1  from  a  posterior  probability  of  0.924 
to  0.106  when  using  the  parameter  prior  distribution.  However,  using 
the  parameter  posterior  distribution  gave  consistent  results  as  shown 

i n  Tab  1 e  4.3* 


Thus,  the  model  discrimination  procedure  is  very  sensitive 
to  the  parameter  prior  when  using  the  parameter  prior  distribution  to 
compute  the  expected  likelihood.  In  fact,  as  can  be  seen  from  relation 
2.13,  the  expected  likelihood  of  the  model  is  approximately  propor- 
tional  to  |C^|  for  large  C^.  Hence,  if  the  prior  for  a  particular 

model  is  made  somewhat  more  vague,  keeping  the  same  priors  for  the 
remaining  models,  the  posterior  probabilities  will  be  drastically 
affected.  One  might  obtain  almost  any  set  of  results  by  specifying 
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a  somewhat  different  vague  prior.  This  result  lends  unequivocal  sup¬ 
port  for  using  the  posterior  rather  than  the  prior  parameter  distri¬ 
bution  in  computing  the  expected  likelihoods. 

Hitherto,  the  standard  Bayesian  model  discrimination 
procedure  has  been  to  employ  parameter  prior  distribution  to  describe 
the  parameter  space  for  the  likelihood  function  [Reilly  (1970),  Box 
and  Hill  (1967),  Hill  (1966),  etc].  It  might  therefore  seem  surpris¬ 
ing  that  no  one  has  experienced  these  problems  in  the  face  of  evident 
inconsistency.  The  answer  probably  is  that  no  one  has  tried  to 
incorporate  these  vague  priors  in  the  model  discrimination  problems. 
The  usual  practice  has  been  to  get  around  the  problem  by  using  part 
of  the  data  to  obtain  a  prior  distribution  of  the  parameters  and  then 
continue  further  by  using  this  prior. 

Box  and  Hill  (1967)  considered  the  following  four  polynomial 
models  for  sequential  model  discrimination. 


Model  1 . 

Model  2. 
Model  3* 
Model  4. 


y,  =  a,x 


y2  =  a2  +  b2x 


y3  =  a3  +  b3x  +  c3x 


Vi,  *  >v  +  V 


The  data  was  generated  from  model  3  with  known  error  variance  a  =  1 
and  a^  =  b^  =  c^  =  1 .  To  start  with,  all  the  models  were  considered 
equi-likely.  Table  4.4  gives  the  results  due  to  Box  and  Hill  [Hill 
(1967)]  and  Table  4.5  gives  the  results  using  the  same  data  but 
employing  posterior  distribution  of  parameters. 
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It  is  not  clear  as  to  what  procedure  Box  and  Hill  used  in 
arriving  at  the  posterior  probabilities  for  the  initial  five  data 
points.  It  is,  however,  evident  that  they  could  not  have  started 
with  an  improper  prior  because  of  the  problems  associated  with  it 
while  using  the  parameter  prior  distribution  in  computing  the 
expected  likelihoods.  Examining  the  posterior  probabilities  after 
the  initial  five  data  points  from  Tables  4.4  and  4.5,  it  appears 
that  Box  and  Hill  used  the  first  five  data  points  to  generate  a 
prior  and  subsequently,  they  updated  the  model  probabilities  by 
using  the  prior  distribution  of  the  parameters  for  every  additional 
data  point  generated.  If,  this  were  the  case,  all  the  subsequent 
model  discrimination  is  really  being  carried  among  the  specific  members 
of  the  general  models  obtained  after  incorporating  the  first  five 
data  points.  The  use  of  posterior  distribution  of  parameters  also 
implies  d i scr imi nat i on  among  specific  members  of  the  families,  but 
these  members  represent  the  "best"  member  of  each  family  at  a  given 
point  in  the  data  gather i ng  process .  However,  the  results  of  model 
discrimination  for  the  given  problem  are  not  significantly  different 
in  Table  4.5  by  using  the  posterior  distribution  of  parameters  from 
those  reported  by  Box  and  Hill  in  Table  4.4. 

Problems  associated  with  the  use  of  parameter  prior  distri¬ 
bution  become  still  more  severe  for  the  non-linear  models  because  of 
the  additional  approximations  due  to  linearization.  Reilly  (1970) 
suggests  the  use  of  posterior  parameter  means  but  the  prior  variance- 
covariance  matrix  for  obtaining  expected  liklihoods.  This,  however, 
still  leaves  the  problem  of  hypersensitivity  of  the  results  to  the 
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Table  4.4 

Results  When  Discriminating  Among  the  Four  Polynomial  Models 
Using  Prior  Parameter  Distribution 

J- 

In  Computing  the  Likelihood 


n 

X 

y 

P1  ,  n 

P2,n 

P3,n 

P4,n 

0 

.2500 

.2500 

.2500 

.2500 

1 

0.0 

1 .431 

2 

1  .0 

2.576 

3 

2.0 

7.540 

4 

3.0 

11.765 

5 

4.0 

20.442 

.0024 

.0058 

.6583 

.3335 

6 

0.0 

1  .837 

.0009 

.0012 

.8777 

.  1202 

7 

0.0 

0.140 

.0017 

.0022 

.7471 

.2490 

8 

0.0 

1  .686 

.0007 

.0013 

.9038 

.0942 

9 

0.0 

1  .714 

.0002 

.0009 

.9707 

.0282 

[Hill  (1966)]. 
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Table  4.5 

Results  When  Discriminating  Among  the  Four  Polynomial  Models 
Using  Posterior  Parameter  Distribution 


X 

In  Computing  the  Likelihood 

y  Pl,n  P2,n 

P3,n 

P4,n 

.2500 

.2500 

.2500 

.2500 

0.0 

1  .431 

1  .0 

2.576 

2.0 

7.540 

3.0 

11.765 

4.0 

20.442 

.0006 

.0006 

.6652 

.3336 

0.0 

1.837 

.0002 

.0001 

.9099 

.0898 

0.0 

0.140 

.0003 

.0002 

.8342 

.1653 

0.0 

1 .686 

.0001 

.0001 

.9486 

.0512 

0.0 

1  .714 

.0001 

.9866 

.0133 

9 
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prior  variance-covariance  matrix  for  vague  or  near  vague  priors. 

One  disadvantage  of  using  parameter  posterior  distribution 
is  the  inability  to  carry  out  the  sequential  discrimination  by  using 
intermediate  parameter  distributions  as  the  priors  for  the  next  step 
and  update  the  model  posterior  probabilities  by  using  few  observations 
at  a  time.  Instead,  the  prior  distribution  of  the  parameters  before 
any  data  are  taken,  is  used  and  all  the  data  has  to  be  used  together 
for  updating  the  model  probabilities.  This  could  involve  somewhat 
more  computations  than  if  parameter  prior  distribution  procedures  were 
used  in  sequential  model  discrimination.  However,  with  the  more 
efficient  expressions  developed  here  for  the  expected  likelihood  in 
terms  of  the  sufficient  statistics  of  the  data,  the  computations 
involved  for  one  data  point  at  a  time  or  all  the  data  points  together 
are  much  the  same.  In  fact,  with  the  improper  uniform  prior,  as  often 
is  the  case,  computations  using  parameter  posterior  distribution  are 
almost  trivial  compared  with  those  required  in  sequential  updating 
using  the  prior  distribution  for  the  parameters, 

4 . 2  Rationalized  Bayesian  Model  Discrimination  Procedure 

The  recommended  model  discrimination  procedure  is  outlined 

below: 

a)  Assume  the  errors  to  be  normally  distributed  and  proceed  with  the 
discrimination  process.  Examine  the  actual  errors  in  the  predictions 
with  the  best  model.  Existence  of  any  correlation  of  the  errors  with  the 
independent  variables  over  the  entire  operating  space  indicates  the 
ineptness  of  the  model  and  that  it  needs  to  be  more  sophisticated. 
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b)  Ordinarily,  the  discrimination  analysis  may  be  carried  out  assuming 
the  error  variance  to  be  known  and  estimating  the  error  variance  from 
the  actual  sum  of  errors  squared.  For  better  discrimination,  it  is 
desirable  to  use  the  updated  unbiased  estimates  of  the  error  variance 
in  computing  the  parameter  posterior  distributions  and  the  expected 
likelihoods.  However,  for  more  appropriate  results,  the  error 
variance  should  be  treated  as  unknown  and  its  joint  distribution  with 
the  parameters  should  be  obtained.  This  joint  distribution  is  then 
used  to  compute  the  expected  likelihoods  instead  of  using  the  joint 
distribution  of  only  the  parameters  [see  Lindley  (1965),  Raiffa  and 
Sch 1  a i fer  (1961)]. 

c)  It  is  desirable  to  use  the  parameter  posterior  distribution  of 
parameters  to  obtain  consistent  results  even  when  vague  or  nearly 
vague  priors  are  used.  Furthermore,  the  use  of  the  posterior  dis¬ 
tribution  of  parameters  avoids  oversensitivity  of  results  to  the 
priors  for  non-linear  models  also. 

d)  To  obtain  consistent  results  while  using  the  posterior  distribution  of 
parameters  in  the  sequential  discrimination,  the  sufficient  statistics 

of  the  entire  data  must  be  employed  in  computing  the  expected  likeli¬ 
hood.  The  model  discrimination  should  then  be  carried  out  with  the 
original  priors  and  this  expected  likelihood. 
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5.  APPLICATION  EXAMPLES 


5 .1  Petroleum  Reservoir  Estimation  from  Material  Balance 

This  chapter  is  intended  to  be  an  exercise  and  test  for  the 
techniques  of  parameter  estimation.  Estimation  techniques  for  the 
parameters  of  a  petroleum  reservoir  material  balance  equation  (alge¬ 
braic)  and  for  the  parameters  in  a  set  of  non-linear  chemical  kinetic 
differential  equations  are  considered  here.  We  have  assumed  that  the 
model  form  is  known  for  certain  and  that  the  only  uncertainty  lies 
in  the  values  of  the  parameters. 

The  nomenclature  employed  for  the  material  balance  equation 
for  a  petroleum  reservoir  is  the  one  recommended  by  the  Society  of 
Petroleum  Engineers  in  1956  and  is  detailed  by  Amyx,  Bass  and  Whiting 
(I960).  The  material  balance  at  any  time  t,  for  a  general  case  of 
a  petroleum  reservoi r  with  an  original  gas  cap  and  water  drive  yields 


N  [B 
P 


t 


+  B  (R 

g  p 


=  N  [  ( B 


t 


R  .  )  +  W  -  W.  -  G. B  . ] 
si  p  1  1  g 1 


B 


-  B  .)  + 


t  i 


.  .  c  (C,  +  S  C  )AP 
1 1  1  -S  f  w  w 


w 


mB 


(B  -  B  )]  +  W 
B  .  g  g  1  e 


(5.1) 


where 

B  -  gas  formation  volume  factor  at  time  t 
9 

B  .  -  qas  formation  volume  factor  at  t  =  0 

gi 

B^  -  total  oil  formation  volume  factor 

B  .  -  initial  total  oil  formation  volume  factor 
1 1 

C^-  -  rock  formation  compressibility 
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-  water  compressibility 
G.  -  cumulative  gas  injected 

m  -  ratio  of  initial  gas-cap-gas-reservoir  volume  to  the 

initial  oil  in  pi  ace 

N  -  initial  oil  in  pi  ace 

Np  -  cumulative  oil  produced 

AP 1  -  reservoir  pressure  decline  since  t  =  0 

Rp  -  producing  gas-oil  ratio 

R  .  -  initial  so  1 ut i on-gas-oi 1  ratio 
s  i  3 

S  -  water  saturation 
w 

W  -  cumulative  water  influx 
e 

W,  -  cumulative  water  injected 
Wp  -  cumulative  water  produced 

The  term  on  the  left  of  equation  (5.1)  represents  net  produc¬ 
tion  from  the  reservoir  in  barrels  at  the  reservoir  conditions  and 
shall  be  denoted  by  F.  An  excellent  treatment  on  the  evaluation  of 
the  water  encroachment  term  is  given  by  Craft  and  Hawkins  (1959). 

Van  Everdingen  and  Hurst  (19^9)  developed  a  procedure  to  evaluate 
by  superposition  of  the  contributions  due  to  the  pressure  declines 
between  all  the  successive  time  intervals  which  may  be  expressed  as 

We  ■  C,  l  *P-a<V  re/rJ 

where  the  summation  extends  over  all  the  time  intervals  until  time  to 
AP  is  the  average  rate  of  pressure  decline  at  a  particular  time  point 
in  the  summation  and  is  usually  computed  by  taking  the  average  of  the 
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pressure  declines  in  the  preceding  and  the  following  time  intervals. 

Van  Everdingen  and  Hurst  (19^9)  and  Craft  and  Hawkins  (1959)  have 

tabulated  Q.  as  a  function  of  the  ratio  of  aquifer-reservoir  radii 

r  /r  and  the  dimensionless  time  t-  =  (At^.t1  where  t1  is  the  actual 
e  w  D  D 

time  in  intervals  starting  from  the  interval  under  consideration  and 
up  to  the  total  time  t.  At^  is  a  parameter  to  convert  the  actual 
time  to  the  dimensionless  units  and  is  simply  a  constant. 

Equation  (5.1)  may  be  simplified  by  neglecting  the  rock  and 
water  expansion  as  for  saturated  reservoirs  and  becomes 

B  . 

F  =  NE  +  Nm  — E  +  C.  Y  AP.Q  (t  ,  r  /r  )  (5.2) 

o  B.g  1  L  D  ew 

g  i 

Here,  E  =  B  -  B  .  refers  to  the  oil  expansion  and  E  =  B  -  B  .  is 
o  t  1 1  r  g  g  g  i 

the  gas  expansion  in  bringing  them  from  the  reservoir  conditions  to 

the  reference  conditions.  F,  E  ,  E  ,  B  .  ,  B  .  and  AP  in  (5*2)  can  be 

ogtigi 

obtained  from  the  actual  field  measurements.  Q(t-,  r  /r  )  may  be  read 

D  e  w 

from  the  tables  provided  At-,  and  r  /r  are  known.  Thus  we  have  five 

r  Dew 

unknowns  in  (5.2),  viz.  the  initial  oil  in  place  N,  the  ratio  of  initial 

gas  cap  volume  to  the  initial  oil  rock  volume  m,  C^,  At^  and  the  ratio 

of  aquifer-reservoir  radii  r  /r^.  Semi -emp i r i ca 1  methods  are  available 

for  the  estimation  of  N,  m  and  AtQ  [Amyx ,  Bass  and  Whiting  (i960)]  but 

are  seldom  reliable.  It  is  therefore  desirable  to  estimate  these 

parameters  from  the  production  history  of  the  reservoir.  As  can  be 

seen  N,  Nm  and  C..  are  the  linear  parameters  and  At-  and  r  /r  are  the 
’I  Dew 


non-linear  parameters. 
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Since  the  precision  of  the  field  data  is  usually  poor, 
estimation  of  all  the  five  parameters  using  the  general  case  of 
equation  (5*2)  may  lead  to  rather  unrepresentative  estimates  of  the 
parameters.  Besides,  it  involves  more  computations  particularly  in 
the  non-linear  search.  Hence,  any  information  that  simplifies  the 
general  equation  should  be  incorporated  to  reduce  the  model  to 
fewer  parameters.  Now,  we  consider  the  various  reservoir  situations 
in  the  order  of  complexity.  However,  it  is  implicit  that  all  the 
following  analysis  neglects  the  rock  and  water  expansion. 


5.1.1  No  Water  Drive,  No  or  Known  Original  Gas  Cap 

When  no  original  gas  cap  is  present  (5-2)  yields 


(5.3) 


This,  a  trivial  case,  and  parameter  estimation  constitutes 
obtaining  the  least  squares  estimate  of  N.  From  a  record  of  n  time 
intervals,  we  obtain  an  nxl  design  matrix,  such  that 


X  =  (x. .}  =  (E  }. 

—  1  j  01 


Compute 


A 

A 

N 


xTx 

A-'xTf 


S  =  (f  -  Ne  )  1  (f  -  Ne  ) 

—0  —  — o 

where  f  is  the  vector  containing  cumulative  production  F  from  the 


reservoir  for  the  n  time  intervals  and  e  is  the  vector  containing 


-o 


E  . 
o 


The  unbiased  variance  of  the  model  is  estimated  by 
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2 

a 


m=  1 


The  variance  of  the  estimate  of  oil  in  place  N  is  then  the 
only  element  of  the  variance-covariance  matrix  C  defined  by 


C 


a2  A 


-1 


a 


(N) 


C  [  1  ,  1] 


95%  confidence  interval  of  N  =  +_  2a  ^  . 

The  case  of  known  gas  cap  is  identical  to  the  above  except 

B  . 

that  E  is  replaced  by  E^.  =  E  +  m  - — —  E  . 
o  r  1  t  o  Bog 

9> 

5.1.2  No  Water  Drive,  Unknown  m  and  N 


F  =  NU  +  m  E  ) 

0  Bgi  9 

or 

F  E 

—  =  N  +  G  ^  (5.4) 

E  o 

o 

Havlena  and  Odeh  (1963)  have  suggested  the  following  two 
methods  of  estimating  m  and  N 

Bti 

(a)  Assume  an  m  and  plot  Eq  +  m  - —  E  =  E  vs  F.  Adjust  m  so  that 

g«  9 

a  straight  line  through  the  origin  is  obtained, 

(b)  Plot  F/E  vs.  E  /E  so  as  to  obtain  a  straight  line  with  N  as  an 

o  g  o  3 

intercept  and  G  as  a  slope. 

To  employ  statistical  methods  of  modeling  and  parameter  esti¬ 
mation,  the  following  approach  is  recommended  for  the  above  methods. 
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(a)  Consider  the  non-linear  model 


F  =  N(E  +  m  E  )  +  C'  (E  +  m 
o  B  .  g  o 

9  i 


(5-5) 


where  C1  is  a  nuisance  parameter  and  physically  represents  the  curv¬ 


ature  of  F  vs.  E  plot.  Non-linear  estimation  is  employed  to  estimate 


m  using  any  of  the  techniques  prevalent  in  literature  [Marquardt  (1953), 
Powe  1  1  (1964 )  ]  . 


However,  the  criterion  for  the  parameter  estimation  is  to 


minimize  the  sum  of  errors  squared  subject  to  keeping  the  nuisance 
curvature  parameter  C1  to  an  acceptably  small  value.  For  each  value 
of  m,  N  and  C1  are  estimated  by  linear  least  squares.  The  final 
estimates  of  N  are  obtained  by  the  linear  least  squares  using  C1  =  0 
and  the  estimated  value  of  m.  The  variance  of  the  fit  and  the  con¬ 
fidence  intervals  of  the  linear  parameters  may  be  estimated  as  in 
section  5.1.1.  However,  the  variance-covariance  matrix  of  the  linear 
parameters  shall  now  include  all  the  uncertainties  associated  with  the 
non-linear  parameter  m. 

(b)  This  is  a  case  of  a  simple  linear  model. 


Let 


F 


o 


form  X  such  that 


i  =  1  ,  n 


E 


9 


(5.6) 


o 
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The  parameters  and  the  variance-covariance  matrix  is  obtained 
as  in  section  5.1.1.  The  variance  of  N  and  G  is  given  by 


2 

°(G) 


=  C  [  1  ,  1  ] 

=  £[2,  2] 


5.1*3  Water  Driven  Reservoirs,  Two  Unknowns 

With  water  drive  present  but  no  or  known  original  gas  cap, 
the  MBE  reduces  to 


where 


F 


N  +  C 


I  AP.Q(t0, 


r  /r  ) 
e  w 


E  + 

o 


m 


fti. 

B  . 

g  > 


E 

g 


(5.7) 


There  are  four  parameters  in  the  above  model,  two  linear,  N 

r 

0 

and  C  and  two  non-linear,  — and  At~. 

r  D 

w 

Application  of  analytical  techniques  for  the  estimation  of 
non-linear  parameters  (r  /r^)  ar|d  Atp  is  limited  by  computations  of 
Q(t^,  r  /r^),  t^ie  rigorous  calculations  of  which  require  vast  computa¬ 
tional  efforts  at  each  step.  It  is  therefore  desirable  to  approximate 
it  by  an  empirical  polynomial  from  the  tabulated  values.  One  such 

computer  program  'WDR1  that  evaluates  Q(tp,  r  /r  )  in  the  range  of 
r  r 

—  =  5  "  10  with  an  accuracy  of  +_  2%  and  for  —  =  1.5  -  5,  within  +_  5% 
w  w 

is  given  in  the  appendix. 
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The  statistical  parameter  estimation  procedure,  then,  consists 

r 

of  employing  a  two-dimensional  non-linear  search  for  — and  At-  and 

r  D 

w 

simultaneously  estimating  the  linear  parameters  N  and  C  as  indicated 
by  the  following  steps: 

(a)  At  a  particular  iteration,  choose  or  arrive  at  a  value  of  r  /r  and 

e  w 

AtD* 

(b)  Compute  the  linear  parameters  N  and  C  for  the  known  values  of 
rg/rw  and  AtD  by  the  linear  regression  described  in  section  5 °  1  .  1  . 

(c)  Compute  the  sum  of  errors  squared  which  now  includes  the  errors 
due  to  the  linear  as  well  as  the  non-linear  parameters  and  of  course 
the  measurement  errors. 

(d)  Choose  the  set  of  r  / r  and  At-,  that  leads  to  minimum  sum  of 

e  w  D 

errors  squared  fit  for  the  resulting  linear  model. 

Experience,  however,  has  shown  that  minimizing  the  sum  of 
errors  squared  alone  is  not  adequate.  It  is  also  required  that  the 
errors  be  uncorrelated  over  the  time  space.  To  avoid  obtaining  a 
misleading  fit,  it  is  proposed  to  use  the  following  five-parameter 
mode  1 


F 

—  M  -i-  r 

y  AP.Q(tn,  r  /r  ) 

L  Dew 

+  S 

I  AP.Q(tD>  re/rw) 

E  N  C1 

E 

E 

t 

t 

t 

(5.8) 

The  criterion  now  is  to  obtain  r  /r  and  At-  that  minimizes 

e  w  D 

the  sum  of  errors  squared  subject  to  keeping  C^  within  a  small  accep- 


tab 1 e  value  e  ,  i  .e . 
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min.  S 


C2I  <  e 


The  variance-covariance  matrix  for  the  linear  parameters 
now  includes  the  uncertainties  associated  with  the  non-linear 
parameters . 


5-1. 4  General  Case,  Unknown  m,  N  and  Water  Drive 

Havlena  and  Odeh  (1964)  reported  maximum  difficulties  in 

this  case  owing  chiefly  to  the  necessity  of  obtaining  derivatives  of 
Bti 

- —  E  for  use  in  the  existing  straight  line  method.  This  necessitates 
9'  9 

the  availability  of  extremely  accurate  data  for  any  reliable  results 
to  be  obtained. 

However,  the  non-linear  estimation  of  parameters  described 
here  does  not  require  any  such  derivatives  and  is  expected  to  perform 
better  than  the  straight  line  method  proposed  by  Havlena  and  Odeh 
(1963).  The  generalized  MBE  is 


-= N  +  G  if  + 

o  o 


C'l  AP.Q(AtD) 


(5.9) 


o 


This  is  a  five-parameter  model  with  N,  G  and  C1  as  the  linear 
parameters  and  r^/r^  and  At^  as  the  non-linear  parameters.  The  para¬ 
meter  estimation  and  their  variance-covariance  matrix  may  be  obtained 
as  in  Section  5.1*3,  except  that  now  the  design  matrix  formed  will  be 
of  the  order  nx3  for  the  linear  parameters.  Alternatively,  a  six- 
parameter  model  as  in  Section  5.1.3,  may  be  employed. 


7b 


F  E 

_ -  M  +  n  9  +  r 

I  AP.Q(AtD) 

4-  r 

l  AP . Q. ( A t  D ) 

t  t  1 

0  0 

E 

0 

C2 

E 

0 

The  parameters  are  estimated  by  minimizing  the  sum  of  errors  squared 
subject  to  keeping  within  a  small  acceptable  limit. 

5.1.5  Incorporating  Prior  Information 

Frequently,  oilmen  have  some  prior  knowledge  of  the  oil  in 
place,  size  of  the  gas  cap  and  reservoir  characteristics  from  other 
exploration  activities  such  as  volumetric  estimates  of  the  oil  in  place 
or  geologists'  personal  estimates  based  on  their  past  experience, 
Bayesian  statistical  methods  for  parameter  estimation  enable  us  to 
combine  the  prior  information  with  the  new  production  data  in  arriv- 
i ng  at  the  final  results. 

If  the  prior  beliefs  are  only  vague  in  the  sense  that  no 
particular  value  for  the  oil  in  place  seems  more  likely  than  the 
others,  infinite  variance  may  be  assigned  to  the  various  parameters 
and  such  a  prior  knowledge  is  called  the  uniform  improper  prior.  Theo¬ 
retically,  such  an  i nformat i on  1  ess  situation  can  never  be  encountered 
in  oil  exploration,  since  anyone,  for  example,  can  assert  that  the 
oil  in  place  cannot  be  negative.  Nevertheless,  the  uniform  improper 
prior  can  be  employed  if  the  prior  beliefs  are  rather  vague. 

The  procedure  of  the  parameter  estimation  and  determining 
the  confidence  intervals  of  the  parameters  for  the  general  case  now 
becomes  as  follows: 

(a)  Assign  the  expected  values  to  the  various  linear  parameters  and 

call  it  b  . 

T> 
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(b)  Next,  is  the  more  difficult  task  of  quantifying  the  strength  of 
one's  beliefs  or  estimates.  In  this  case,  a  3  x  3  variance-covariance 

C 1  2  C  1  3 
C22  C23 
C32  C33 

where  C  ^ and  are  respectively  the  marginal  variances  of  N, 

G  and  C^.  C.^  =  C^,  C ^ 3  =  C.^>  etc’  are  the  covariances 

or  a  measure  of  correlation  between  the  two  parameters.  Usually,  it 
is  hard  to  assign  any  correlation  between  the  various  parameters  on 
physical  grounds.  Hence,  it  is  usual  to  make  all  the  off-diagonal 
elements  of  to  be  zero  implying  a  prior  belief  of  no  correlation 
between  the  various  parameters. 

The  variances  C^,  may  be  assigned  with  relative 

ease.  For  example,  a  geologist  may  claim  that  the  volumetric  estimates 
of  oil  in  place  are  23.2  MMB  and  he  feels  about  95%  sure  that  the 
actual  oil  in  place  will  not  be  off  from  this  figure  by  more  than  +4 
MI1B,  Assuming  normal  distribution  of  his  beliefs,  as  we  have  done 
for  all  the  measurement  errors  and  the  parameters  in  our  statistical 
analysis,  4  MMB  corresponds  to  two  standard  deviations  and  thus  prior 

=  2,  giving  =  a  ^  =  4.  Of  course  making  infinite  or  very 

large  values  implies  that  little  is  known  about  the  parameter  Cp 


matrix  has  to  be  specified. 
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(c)  Estimate  or  assign  a  value  to  the  error  variance  arising  chiefly 

2 

due  to  the  measurement  errors  and  call  it  , 

(d)  Form  the  n  x  3  design  matrix  X_  from  the  production  data  at  n 
time  intervals  as  described  in  earlier  cases. 


Define  A  = 

L  = 

A  = 

-o 

where  y  is  a  vector  containing  n 
of  the  model,  i.e. 

o 

(e)  Estimate  the  updated  values 
covar i ance  matrix 


X 


observed  values  of  the  left-hand  side 


of  the  parameters  and  variance- 


A  =  A  +  A 

— n  — o  — 


b  =  A  ^  {A  b  +  A  b) 

— n  — n  —o—o - 


c  = 


S 

n-m 


m=3 


A  revised  estimate  of  error  variance  may  be  obtained  by  computing  a 

2  2. 

weighted  average  of  a  and  according  to  certain  beliefs  in  the 

2 

prior  estimate  o  in  terms  of  equivalent  observations.  Call  this 


revised  estimate  <3. 


Variance-covariance  matrix  is  now 


C  =  a  2A_1 
— n  F  — n 


Variance  of  N,  G  and  C  are  the  respective  diagonal  elements  of  C 


— n 


77 


This  approach  can  be  used  to  update  the  estimates  of  N  and  G 
as  and  when  new  f i el d  data  is  obtained  without  resorting  to  the  treat¬ 
ment  of  entire  earlier  data  which  has  already  been  processed.  The 
parameter  estimates  and  the  variance-covariance  matrix  of  the  previous 
data  is  treated  as  the  prior  for  the  new  analysis. 

5 . 2  Petroleum  Reservoir  Es t i ma t i on~~F i el d  Examples 

Examp  1 e  5.1 

This  is  a  case  of  water  drive  without  any  gas  cap.  The 
example  is  taken  from  Van  Everdingen  et  al  (1953) •  An  independent 
volumetric  estimate  of  oi 1  in  place  is  24-25  MM  EL 

Since  in  equation  (5.8)  was  found  to  be  very  small  for  the 
case  of  infinite  water  drive,  the  model  employed  was 

F 

eT 

o 

It  was  decided  to  estimate  the  parameters  by  using  partial 

data  and  compare  the  results  as  more  and  more  data  was  employed. 

Parameters  were  estimated  by  minimizing  the  sum  of  errors  squared. 

Since  re/rw  =  00  was  found  to  give  markedly  superior  fit  as  compared 

to  some  of  the  other  sets  of  r  /r  tried,  one-dimensional  non-linear 

e  w 

search  was  employed  for  At^  with  r  / rw  =  The  following  polynomial 

fit  was  obtained  from  the  tabulated  values  of  Q.  ( t  D ,  r^/r^)  for 

r  /r  =  oo  [Craft  and  Hawkins  (1959)]  and  was  used  here, 
e  w 


l  AP.Q(tD,  re/rw) 
E 
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Q ( t D )  =  -  .251233  +  1  .50397^tD1/2  +  .2955831 tp 

9  -7  7  —in  Zi 

-  , 000 1 426^3 1 D  +  1,512  X  10  't^  -  1.1496  x  10  tQ 

This  polynomial  was  found  to  give  Q(t^)  correct  within  I .5%  over  the 
range  of  tD  =  10  to  700. 

The  reservoir  started  production  in  1942  and  one  data  point 

was  available  marking  the  end  of  each  quarter  up  to  the  last  quarter 

of  1950.  Partial  data  sets  were  considered  until  the  end  of  1944, 

1946,  1948  and  the  entire  data  until  the  end  of  1950,  Table  5.1  gives 

the  estimates  of  the  initial  oil  in  place  N  and  the  variance  of  these 
2 

estimates  a ^  for  the  various  sets  of  partial  data.  The  estimates 
of  the  original  oil  in  place  are  reasonably  good  even  when  only  the 
first  three  years  of  data  have  been  used.  At^  =  14  was  found  to 
give  the  minimum  sum  of  errors  squared  using  the  entire  data  set  the 
same  as  reported  by  Van  Everdingen  et  al,  after  a  complex  analysis 
of  the  variations  in  N,  calculated  at  each  time  interval  from  the 
material  balance  equation.  The  original  oil  in  place  was  found  to  be 
25.92  HUB  and  £  =  0,000193.  Variance-covariance  matrix  for  N  and 


£  = 

6.57  x 

ID'2 

-1.38  X 

10 

-I.38  X 

O 

1 

3,24  x 

10 

This  variance-covariance  matrix,  however,  also  includes  the  uncertainties 

associated  with  non-linear  parameters  At~  and  r  /r  . 

D  e  w 

2  2 

The  variance  of  fit  is  a  =  0,23  MM B  ,  thus  the  standard 
deviation  of  the  estimates  of  oil  in  place  0  ^  =  /C [ 1 , 1 ]  =  /0 . 0657 
=  0.256  MMB  and  95%  confidence  interval  for  N  =  +_2a,^ 


=  +  0.513  MMB. 
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Table  5.2  gives  the  detailed  calculation  using  the  entire  data  set  and 
Atp  =  14.  The  average  pressure  decline  AP ,  is  computed  by  averaging 
the  pressure  decline  in  the  preceding  and  the  following  intervals. 

The  dimensionless  time  tp  =  ( At  D ) » t '  is  obtained  by  computing  the 
actual  time  t1  in  intervals  starting  from  the  given  interval  and 
up  to  the  final  interval  of  the  data.  The  polynomial  given  above  was 
used  to  get  Q(tp) 

Example  5.2 

We  use  this  example  to  illustrate  the  effect  of  prior  infor¬ 
mation  on  the  final  results  of  Example  5.1.  The  volumetric  estimates 
of  the  original  oil  in  place  are  between  24  and  25  MM B  or  a  mean  value 
of  24.5  11MB.  Supposing  now  that  the  estimating  geologist  feels  that 
the  actual  oil  in  place  lies  between  24.5^8  MMBwith  95%  confidence 
or  24.5  12  MMB  with  99%  confidence.  Also,  he  has  no  idea  as  to  the 

value  of  the  linear  water  drive  parameter  ,  to  which  we  arbitrarily 

assign  a  value  of  1,  then 


^  =  [24.5,  1] 


„  _  8  _  12  _  k 
°(N)  "  2  ”  T  " 


0(c)  =  00  @  1 000  (say) 


C  = 

— o 


16  0 


0  10 


Also  assume  that  after  examining  the  measurement  techniques  for  the 
total  cumulative  production  F,  the  petroleum  engineer  assigns 


o E  =  0.5  MMB  in  F/Eq. 


80 


Table  5.1 

Results  for  Example  5.1 

r  /  r  =  °° 
e  w 


Data  Used  Up  To 
and  Including 

Estimates  of 
tD  =  10 

the 

tD 

Original 
=  14 

Oi  1  i 

tD 

n  Place 
=  20 

N 

>)2 

N 

°(N)2 

N 

°(N)2 

\Skk 

26.8 

0.59 

26.1 

0.31 

25.5 

0.51 

1946 

24.4 

0.50 

25.3 

0.27 

26.1 

0.35 

1948 

23.9 

0.22 

25.1 

0.13 

26.2 

0.16 

1950 

24.7 

0.09 

25.9 

0.07 

27.  1 

0.08 
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Table  5.2 

Details  of  Calculations  for  Example  5  - 1 


r  /r  =  °°,  At_  =  14 
e  w  ’  D 


Year 

Quarter 

Boundary 
Pressure 
@  8075  ft. 

Avg .  AP 
for 

Quarter 

- - - 

£APQ(tD)  x  10 

E 

0 

F  in"6 

—  x  10 

0 

T 

o 

0 

3793 

1942 

1 

3788 

2.5 

2 

3774 

9.5 

3 

3748 

20.0 

4 

3709 

32.5 

1943 

1 

3680 

34.0 

0.3813 

.004225 

99.1727 

2 

3643 

33.0 

0.2493 

.010032 

74.3878 

3 

3595 

42.5 

0.2078 

.01759 

65.8676 

4 

3547 

48.0 

0.2093 

,02432 

66.0182 

1944 

1 

3518 

38.5 

0.2267 

.02951 

69.6051 

2 

3485 

31  .0 

0.2443 

.03438 

73.2461 

3 

3437 

40.5 

0.2583 

.04002 

75.5010 

4 

3416 

34.5 

0,2736 

.04534 

78.7561 

1945 

1 

3379 

29.0 

0.2786 

.05233 

78.7521 

2 

3358 

29.0 

0.2989 

.05644 

83.0141 

3 

3338 

20.5 

0.3241 

.05922 

88.7914 

4 

3329 

14.5 

0,3473 

.06204 

92.3293 

1946 

1 

3324 

7.0 

0.3754 

.06347 

98.0065 

2 

3319 

5.0 

0.4042 

.06455 

103.9530 

3 

3302 

1 1  .0 

0.4155 

.06835 

106.5307 

4 

3292 

13.5 

0.4359 

.07056 

110.9557 

1947 

1 

3276 

13.0 

0.4509 

.07354 

1 13.8909 

2 

3243 

24.5 

0.4399 

.08117 

111.5379 

3 

3206 

35.0 

0.4359 

.08827 

1 10.2903 

4 

3184 

29.5 

0.4462 

.09272 

112.3273 

1948 

1 

3165 

20.5 

0.4517 

.09809 

113.4767 

2 

3135 

24.5 

0,4506 

.10508 

1  12.9093 

3 

3108 

28.5 

0,4546 

.11118 

113.7178 

4 

3095 

20,0 

0,4706 

.11429 

116.9875 

1949 

1 

3086 

1  1  .0 

0.4863 

. 1 1721 

119.4385 

2 

3103 

-  4.0 

0.5524 

.10877 

132.8372 

3 

3125 

-19.5 

0.6064 

.10380 

143.2210 

4 

3120 

-  8.5 

0.6322 

.10401 

148.1808 

1950 

1 

3H5 

5.0 

0.6582 

.10423 

152,7882 

2 

3114 

3.0 

0.6842 

.10444 

157.6309 

3 

3H5 

0.0 

0.7113 

.10444 

162.8010 

4 

3116 

-  1  .0 

0.7380 

.10444 

167.4780 
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Therefore 


Now 


and 


A  .  o  2  C'1 
-o  E  — o 


1.56  x  10 
0 


-2 


a2  c' 


3.2  x  10 
1.36  x  107 


b  =  [  25.926 


A 

— n 


A  +  A 

—  — o 


b  =  A  (A  b  +  A  b) 

— n  — n  —0—0 - 


=  [  25.9197 


C 

— n 


a2  A’’ 

— n 


6.5397  x  10 
-1 .3706  x  10 


-2 


-7 


0 

2.5  x  10 


1 .36  x  107 

6.489  x  1012 
.000193  ] 


.000193  ] 


-1 .3706  x  10 

3.22653  x  10 


Suppose  that  the  geologist  had  somewhat  more  confidence  in  his  volu¬ 
metric  estimates  so  that  his  95%  confidence  intervals  will  correspond 
to  4  MMB,  c  ^  =  2MMB  and  now 

6.248  x  10"2  0 

2.5  x  10"7 


A 

-o 


0 


- 


The  first  element  of  A  is  fourfold  that  of  the  earlier  case,  b  remains 

-o  ’  -o 

the  same.  The  results  are  given  in  Table  5.3  for  various  confidence 
intervals  of  the  prior  estimates  of  oil  in  place. 

Example  5.3 

Frequently,  it  is  not  possible  to  measure  the  boundary 
pressure  and  sometimes,  the  exact  boundary  is  not  even  known.  In 
such  cases,  the  reservoir  average  pressure  is  usually  employed  in  the 
material  balance  calculations  with  certain  expected  reduction  in  the 
accuracy.  It  was  decided  to  examine  the  effect  of  using  average 
reservoir  pressure  in  the  field  case  of  Example  5.1  and  compare  the 
results  for  the  two  cases. 

The  estimation  procedure  is  exactly  the  same  as  in  Example 

5.1  and  the  optimal  results  were  obtained  for  the  radial  drive  with 

t.,  =  31  and  r  /r  =  The  oil  in  place  was  now  23.8  MMB  aqainst  25.9 
Dew 

2 

MMB  obtained  in  Example  5.1  and  the  variance  of  fit  was  a  =  0.2053 

2 

even  lower  than  the  value  of  0.23  MMB  obtained  by  using  the  average 

boundary  pressure.  The  95%  confidence  intervals  were  also  somewhat 

lower  than  those  obtained  in  Example  5*1. 

2 

From  the  values  of  c  and  the  95%  intervals  of  the  original 
oil  in  place,  the  use  of  average  reservoir  pressure  did  not  cause  any 
appreciable  decrease  in  the  accuracy  of  estimation.  Since,  the  volu¬ 
metric  estimates  of  the  oil  in  place  of  24-25  MMB,  are  somewhere  in 
between,  it  is  not  possible  to  assert  the  superiority  of  one  or  the 
other  set  of  values. 

This  indicates  that  the  error  in  estimates  is  more  due  to 
the  errors  in  the  representative  pressure  measurements  than  due  to  the 
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Table  5-3 

Effect  of  Prior  Information  on  the  Final  Estimates 

of  the  0 i 1  i n  P 1  ace 


95%  Confidence  Intervals  Posterior  Estimates  of  Variance 

for  Volumetric  Estimates  Oil  in  Place  of  Final 

Estimates 


00 

25.92600 

.065689 

8.0 

25.91965 

.065397 

4.0 

25.90095 

.064535 

2.0 

25.83083 

.061305 

1  .0 

25.60882 

.051078 

0.5 

25.16505 

.03064 

0.0 

24.50000 

.00000 
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approximation  of  using  average  reservoir  pressure  for  the  boundary 
pressure.  In  fact,  it  is  felt  that  weighted  average  AP ' s  obtained 
from  the  boundary  pressure  and  the  reservoir  pressure  will  lead  to 
more  reliable  results  than  those  obtained  by  using  either  set  of 
pressures  alone. 

5 . 3  Parameter  Estimation  from  a  Set  of  Non-linear 

Differential  Equations 

The  example  chosen  represents  a  contemplated  reaction 
mechanism  for  the  disproportionation  of  freons  on  a  highly  fluorin- 
ated  alumina  catalyst.  Stoi ch i ometr i ca 1 1 y ,  the  reaction  may  be 
described  by  two  reversible  reactions 

2CHC 1 -F  +  CHC  1  F  +  CHC 1 0 

2  2  3  (5.11) 

2CHC 1 F2  +  CHF3  +  CHC 1 2F 

Cavaterra,  Fattore,  and  Giordana  (1967) »  studied  this  reac¬ 
tion  at  200° C .  Their  data  is  given  in  Table  5.*+.  In  interpreting  the 
kinetics  of  this  reaction  system,  Caveterra  et  al  assumed  that  the 
catalytic  disproportionation  could  be  represented  by  the  following 
kinetic  model 

A,  +  uo  a2  +  u, 
k3 
k2 

A1  +  U1  *  \  +  Uo 

k4 

A2  +  Uo  *  A3  +  U1 


(5.12) 
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Table  5.4 

D i spropropor t i onat i on  Reaction  of  Pure  CHCl^F  at  200°C 


Time  Reaction  Components  (Mol.  Fraction) 


(sec) 

CHC 1 2F 

CHC 1 F2 

CHF3 

chci_3 

0.3 

0.850 

0.063 

0.002 

0.085 

0.5 

0.780 

0.106 

0.004 

0.110 

1  .7 

0.435 

0.239 

0.035 

0.291 

2.0 

0.378 

0.255 

0.036 

0.331 

2.3 

0.318 

0.268 

0.047 

0.367 

2.7 

0.289 

0.277 

0.055 

0.379 

3.3 

0.255 

0.286 

0.072 

0,417 

3.7 

0.201 

0.286 

0.078 

0.435 

4.0 

0.176 

0.262 

0.078 

0.484 

6.7 

0.086 

0.215 

0.172 

0.527 

8.0 

0.050 

0.120 

0.235 

0.595 
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where  Aj  =  CHC^F,  A£  =  CHC 1 F2 ,  A  =  CHF  ,  A^  =  CHC 1 3  and 
UQ,  are  the  fluorinating  and  chlorinating  active  centers  of  the 
catalyst . 

The  corresponding  kinetic  model  under  the  steady  state 
conditions  may  be  written  as 


dC. 

dt 

dC; 

dt 

dC. 

6t 

dCj 

dt 


k,c, 


klcl 


k4C2 


k2c, 


k2Cl  +  k3C2 


k3C2  ‘  kltC2 


(5.13) 


Covaterra  et  al,  employed  the  method  of  kinetic  structural 
analysis  developed  by  Wei  and  Prater  (1962)  to  estimate  the  reaction 
rate  constant.  Myint  (1970)  devised  a  smoothing  scheme  and  obtained 
the  reaction  rates  from  the  smoothed  data.  Reaction  rate  constants 
were  then  estimated  from  the  resulting  algebraic  models. 

It  was  decided  to  employ  the  quas i 1 i near i zat i on  procedure  to 
estimate  the  parameters  from  the  system  of  simultaneous  differential 
equations.  The  quas i 1 i near i zat i on  approach  is  a  second-order  iterative 
process  and  once  within  the  convergence  radius,  it  converges  very  fast. 
However,  to  obtain  a  starting  value  good  enough  for  convergence  is  a 
difficult  problem.  An  algorithm  proposed  by  Donnelly  and  Quon  (1970) 
was  employed  to  arrive  within  the  radius  of  convergence  through  the 
technique  of  data  perturbation.  The  essential  idea  is  to  derive  a  new 
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problem  so  that  the  existing  solution  is  within  the  radius  of  conver¬ 
gence  of  this  new  problem  and  the  derived  problem  is  then  perturbed 
towards  the  original  problem  in  a  series  of  steps,  small  enough  to 
ensure  convergence  at  each  successive  step.  Donnelly  (1968)  has 
written  a  FORTRAN  program  for  the  above  algorithm,  which  was  used 
for  the  parameter  estimation  in  this  example. 

With  all  the  four  differential  equations  of  (5*13),  difficul¬ 
ty  was  experienced  in  the  final  stages  of  convergence.  However,  this 
difficulty  disappeared  when  the  last  two  differential  equations  of 
(5.13)  were  integrated  analytically  and  replaced  by  the  algebraic 
equations.  The  results  obtained  are  given  in  Table  5«5  and  compared 
with  the  results  obtained  by  Myint  (1970)  and  Cavaterra  et  al  (1967) 
in  Table  5.6.  The  relative  rate  constant  is  obtained  by  dividing  the 
absolute  rate  constant  by  ky 

The  results  obtained  by  the  quas i 1 i near i zat i on  procedure  are 
reasonably  close  to  those  obtained  by  Cavaterra  et  al  from  the  more 
involved  Wei  and  Prater  ( 1 96 2 )  approach.  Myint  (1970)  also  estimated 
the  order  of  the  reaction  and  therefore  solves  somewhat  different 
problem.  However,  it  is  felt  that  whatever  may  be  the  smoothing 
function,  the  smoothed  data  will  not  retain  the  entire  information 
of  the  raw  data.  Besides,  sometimes  it  might  introduce  a  bias  in  the 
data  since  the  smoothing  function  may  not  track  the  true  behaviour  of 
the  data  generating  process. 

It  might  be  appropriate  to  conclude  this  chapter  with  a  few 
remarks  about  petroleum  reservoir  parameter  estimation.  It  is  desirable 


to  consider  the  most  general  material  balance  model  and  estimate  all 


' 

’ 
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Table  5-5 

Parameter  Estimates  of  Model  (5 » 13) 
Using  Data  of  Table  5.4 


1  .C.  t 

=  0  0|  =  1 

c2  =  o  c3  =  0 

ck  =  0 

Parameter  Est 

i mates 

k 1  =  0,2893 

k2  =  0.2603 

CM 

LA 

O 

II 

OA 

k^  =  0.1461 

Variance-Covariance  Matrix 

0.001 666 

-0.000377 

0.001516 

0.000463 

-0.000377 

0.002595 

0.002408 

-0.003119 

0.001516 

0.002408 

0.006495 

-0,004081 

0.000463 

-0.003H9 

-0.004081 

0.004870 

Comparison  of  the  Parameters  Values  with 
Those  Obtained  by  Myint  (1970)  and  Cavaterra  et  al  (1967) 
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the  parameters  involved.  The  significance  of  each  term  in  the  general 
model  may  then  be  inferred  from  the  marginal  variance  of  the  corres¬ 
ponding  parameter.  If,  however,  the  computations  for  the  general 
case  are  very  large,  and  a  very  strong  a  priori  information  can  be  used 
to  simplify  the  model,  resort  may  be  made  to  the  simpler  cases  dis¬ 
cussed  in  this  chapter. 
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6.  DISCUSSION 


6 . 1  Summary 

The  purpose  of  this  thesis  has  been  to  develop  a  procedure 
for  (i)  discriminating  among  the  rival  models,  and  (ii)  the  design  of 
experiments  to  choose  the  experimental  setting  that  is  expected  to 
result  in  maximum  discrimination  among  the  rival  models. 

Hitherto  the  standard  Bayesian  model  discrimination  proce¬ 
dure  has  been  to  employ  parameter  prior  distribution  to  describe  the 
parameter  space  for  the  likelihood  function.  The  results  of  model 
discrimination  [Chapter  h]  when  using  parameter  prior  distribution 
in  computing  the  expected  likelihood  was  found  to  be  very  sensitive 
to  the  parameter  distribution  for  vague  or  near  vague  priors. 

Constructed  examples  were  used  to  demonstrate  the  relative  inconsistency 
of  using  parameter  prior  distribution  rather  than  the  parameter 
posterior  distribution. 

Also,  the  expected  likelihood  of  the  data  is  simply  the 
average  of  the  likelihood  function  over  the  entire  parameter  space. 

The  posterior  parameter  distribution  being  more  accurate  represen¬ 
tation  of  the  parameter  space,  should  be  used  to  compute  the  expected 
1 i ke 1 i hoods . 

When  errors  are  assumed  to  be  normally  and  independently 
distributed  with  known  error  variance,  the  joint  distribution  of 
the  parameters  is  multivariate  normal.  For  this  case,  new  expressions 
[Chapter  2]  for  the  expected  likelihood  could  be  developed  that  in¬ 
volved  only  the  sufficient  statistics  of  the  data.  The  proposed 
expressions  have  considerable  advantage  over  the  existing  method  in 
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terms  of  the  computational  savings  and  reduced  round  off  errors. 

Frequently,  a  situation  arises  when  more  than  one  model 
can  describe  a  physical  phenomenon  about  equally  well.  In  such 
instances,  it  is  required  to  choose  an  experimental  setting  that 
would  sharpen  the  discrimination  between  the  rival  models.  Model 
discrimination  then  becomes  a  sequential  problem  of  designing  the 
experiments  and  incorporating  the  new  results  for  further  discrimin¬ 
ation.  However,  when  using  the  posterior  parameter  distribution  to 
compute  the  expected  likelihood,  the  entire  data  must  be  used 
against  the  original  parameter  prior  at  each  stage  of  the  sequential 
discrimination. 

A  d i scr imi nat i on  function  R  was  proposed  [equation  (3.21)], 
the  development  of  which  was  based  on  the  concept  of  entropy  as 
a  measure  of  the  uncertainty  associated  with  the  system  as  described 
by  Kulback  (1959)  and  Shannon  (19^8).  The  following  assumptions  were 
made  to  make  the  analysis  feasible:  (i )  the  observations  are 
normally  and  independently  distributed  with  a  known  and  constant 
error  variance;  (ii)  the  model  response  y.  is  linear  in  parameters 
a_.  in  the  neighbourhood  of  the  least  squares  estimates  of  the  parameters 
a_.  ;  (iii)  the  data  is  being  generated  from  among  the  models  under 
cons i derat i on . 

Several  design  criteria  have  been  developed  in  the  past 
few  years;  notable  among  them  being  those  proposed  by  Roth  (1965), 

Box  and  Hill  (1967)  and  Reilly  (1970).  However,  the  criterion 
proposed  in  this  work  evaluates  the  expected  change  in  entropy  or 
the  expected  change  in  the  uncertainty  of  the  situation  and  is  both 


simple  and  exact.  It  also  is  consistent  with  the  attempts  of  earlier 

workers  to  maximize  the  difference  in  the  responses  of  the  various 

models  by  the  inclusion  of  (9n+^,i  "  9n+^)  'n  equation  (3.21)  but 

2  2 

restrains  it  with  the  uncertainty  term  (a  +  an+j). 

Some  illustrative  model  building  examples  have  been  presented 
in  the  field  of  petroleum  reservoir  parameter  estimation  in  Chapter  5. 
Parameter  estimation  from  the  material  balance  equation  has  been 
discussed  for  several  cases  in  order  of  the  complexity.  Simultaneous 
non-linear  search  and  linear  estimation  is  recommended. 

In  addition,  an  exercise  in  parameter  estimation  from  a 
set  of  kinetic  non-linear  differential  equations  has  been  included 
in  Chapter  5«  Qusa i 1 i near i zat i on  procedure  together  with  an 
algorithm  proposed  by  Donnelly  and  Quon  (1970)  was  found  to  give 
results  that  are  in  agreement  with  those  obtained  from  other  methods 
[Cavaterra  et  al  (1967),  My*nt  (1970)]. 

6 . 2  Future  Work 

One  area  of  future  work  would  be  to  develop  expected 
likelihood  expressions  for  the  case  of  unknown  error  variance.  The 
marginal  distribution  for  the  error  variance  is  gamma-2.  The  relative 
intractability  of  gamme-2  is  the  major  source  of  difficulty.  However, 
gamma-2  may  be  approximated  by  a  normal  distribution  with  the  same 
mean  and  variance  as  that  of  the  original  distribution  itself.  The 
resulting  form  is  expected  to  be  more  tractable. 

Expected  likelihood  expressions  involving  only  the  suffi¬ 
cient  statistics  for  a  mu  1 1 i - response  case  would  be  very  desirable. 


* 
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Another  area  for  further  development  is  that  of  the 
sequential  design  of  experiments  for  model  discrimination.  An 
expected  entropy  change  criterion  for  the  case  of  unknown  error 
variance  would  be  very  useful.  The  criterion  needs  also  to  be 
extended  to  the  mu  1 1 i - response  situations  with  the  variance- 
covariance  matrix  for  the  various  responses  both  known  and  unknown. 
A  criterion  to  select  the  experimental  setting  that  takes  into 
account  both  the  model  discrimination  and  the  precise  parameter 
estimation  would  be  ideal. 

There  is  a  great  need  for  a  more  unified  treatment  of 
the  various  techniques  used  in  model  building  and  model  discrimin¬ 
ation.  The  purpose  of  this  treatment  would  be  to  indicate  the 
possible  strategies  and  the  associated  numerical  techniques  that 
one  might  adopt  during  an  investigation.  However,  more  experience 
is  needed  in  the  modeling  techniques  before  such  an  attempt  can  be 
made . 
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APPENDIX 

The  appendix  contains  the  programs  used  in  this  thesis.  The 

notat i on 

used  in  the  following  programs  is  described  here. 

A 

Dispersion  matrix  XtX 

AO 

Prior  estimate  of  the  dispersion  matrix 

B 

Least  squares  estimate  of  the  parameters  (vector) 

BL 

Linear  parameters  (vector)  in  petroleum  reservoir  estimation 

C 

Variance-covariance  matrix 

D 

Determinant,  normalizing  matrix 

DSNM 

A  function  subprogram  to  set  up  the  design  matrix 

ER 

Errors  vector 

FET 

(F  t  ET)  of  sect  ion  5 « 1 

FUCN 

A  function  subprogram  to  compute  the  model  predictions 

1  NX 

Mesh  size  in  the  design  of  experiments 

IT 

Reference  index  for  the  competing  models 

L » LI ,L2 

Expected  likelihoods 

LAM 

A  parameter  in  Marquardt's  algorithm 

M 

Number  of  parameters  in  the  given  model 

MUO 

Prior  estimate  of  the  parameters  (vector) 

N 

Number  of  data  points 

NM 

Total  number  of  competing  models 

PA ,  PW 

Independent  variables 

Q.T 

Q(t~,r  /r  )  of  section  5-1 

D  e  w 

A-2 


R  Negative  of  the  expected  entropy  change 

RD  Aquifer-reservoir  radii  ratio  r  / r 

e  w 

SS  Minimum  sum  of  errors  squared 

SETUP  A  function  subprogram  to  specify  the  prior  parameter 
distribution  for  the  given  model 
T  Time  in  intervals 

TD  Dimensionless  time 

VAR  Error  variance 

X,XX  Design  matrices 

XP  Normal  random  number  generated  from  'SIMULATION' 

Y  Vector  of  observations 

YH  Expected  value  of  the  prediction  from  the  entire  metamodel 


0 


PROGRAM 


'  IN  V' 


a  THIN  PROGRAM  COMPUTES  THE  INVERSE  OF  /I  NON- 
o  SINGULAR  MATRIX  BY  GAUSS -JORDON  PIVOTING. 


VI/7  7[[j]V 

V  Z 4-1 NV  M\I\J\K\P 

[1]  ->3  x  i  (  2  =  p  p;  1  )  A  =  /  pM 

r  2  ]  ->0  ,  pO '  NO  IN  VERSE.  FOUND  ' 

[3]  M*-§  (  1  0  +  p/7)  p  (  ,<$?//)  %J+-1  <P+-\I<-\ \ pM 

[.  4  ]  MlK ,  1  ;  ]«-/'/[  1  ,  K<-K  i  T  I  //[  i  J  ;  1  ]  ;  ] 

[  5  J  Z3«-lct)7,  0pPC;M  >P[  1  ,A'  j 

[  6  ]  //[  ;  1  +  pJ  ]-* — J 

[7]  ->2xii/7"  3  0  >  |  Ml  1  ;  1  ]  t  T  /  |  ,M 

[  3  ]  /P-lcj>l4>[  i  ]  M-  (Jxlll  ;  1  ]  )  °  .  x,7[  1  ;  ]«-//[  1  ;  ]  -:-//[  1  ;  1  ] 

[  3  ]  >4x  ib*I«*I-l 

[10]  ZHllUP] 

V 


PROGRAM  ' DET ' 


o  'DET'  OBTAINS  THE  DETERMINANT  OF  A 

a  SQUARE  MATRIX  BY  OHIO'S  METHOD. 


HDETUY]H 

V  D+DET  A \N \A 

[1J  N<-pA  [  1  ;  J 

[2]  I*-  0 

[3]  D<- 1 

[  4  ]  ->L/15x  i//=  1 

L  5  J  L  A  L  :  I  ■*- 1  +  1 

[  G  ]  D<-D*A  [  1 1 1  ] 

[  7  J  /I  [  J  ;  ]  «-/l  [  I ;  J  */l  [  I  ;J] 

[  8  J 

r  9]  LAD:L<-L  +  1 

[  1  0  1  A  [  iL  1<-A  [  ;  L  ]  -/l  [  I ;  L  \  *A  [  ;  I  ] 

[11]  +LAD*\L<N 

[12]  -+LAL*\I<N-1 

[13]  L/l/?:£«-0x/I [//;//] 

V 


PROG R A M  ' SIMULATION ' 


pi 


ft  'SIMULATION'  GENERATES  NORMAL  RANDOM  NUMBERS 

ft  REFERENCE:  26.2.23,  PC  933,  HANDBOOK  OF  MATH- 

ft  MATICAL  FUNCTIONS ,  (/.S.  DEPARTMENT  OF  COMMERCE. 


MSI  MU  LA TIONtti J V 

V  SIMULA  I ION ; CO ; Cl ; (72 ; 0 1 ; D2 ; D 3 ; P ;P 
L1J  C  0  2  .  615517 

[  2  J  6’ l<-0.  802853 

L  3  J  (72^0.010328 

[  4  J  Jl«-1. 432788 

L  5  J  D 2 +- U  .  18926  9 

[  6  J  D  3^-0.001  308 

[  7  j  P«-((1?1000Q))*10000 

[  8  J  PP+P 

L  9  J  px  iP<Q  .  5 

[  1  o  J  P-*- 1  -P 

[11J  MAT  :'!+(*(  *(P*2)  )  )*0.  5 

[  12  J  XP+-T-  (  ((C0+C1xP)+C2xP*2)t(((1+P1xP)+P2xP*2)+D3xP*3)) 
[  13  j  ->p^px  i?P<0  .  5 
c 14  J  XP+-XP 
[15]  TAT:-+\U 

V 


ft  PROGRAM  '  USJ' 


ft  P//PP  PROGRAM  COMPUTES  THE  EXPECTED  LI KE LI  - 
ft  HOODS  FROM  THE  SUFFICIENT  STATISTICS  OF  THE 
ft  ZMP4  IV//P/7  IMPROPER  UNIFORM  PRIORS  ARE  USED. 


MHSTl UJ V 

V  HST\M\N 
L  1  J  «-  pP 

[  2  J  C/-*-  pXi  ;  1  J 

[3]  Z>(  (  2xo  VAR)*i -Nr  2  )  )x(  0.  5  *  ( M  v  2  ]  )  x  (  *  (  -SS  *  (  2  x  VA  R  )  )  ) 

V 


DROGRAt  1 


'LUMA' 


0 


TRIG  PROGRAM  COMP UTRG  THE  LEAPT  SQUARE  ESTIMATES 
AMD  THE  VAR I A  MCE  -  CO  VAR I A  MCE  MATRIX  FOR  A  LI  SEAR 


MODEL  HUES  ERROR  VAR  I  AM  CE  IS  KSOi/S . 


VLIUM\_{ J  J  V 

V  LI  MM 

Llj  /I  «-(<¥/)+  .  x/ 

r?  ]  h+(tuv  /i)  +  .x(w)  +  .xy 

[3]  A+AiVAR 
[  4  j  C<-I!I  V  A 

L  3  J  SS++/ (Y-X+.*B)*2 

V 


°n 


PROGRAM  MISS' 


'ESS1  COMPUTES  TI!E  EXPECTED  LIKELIHOOD  OF  A  GIVE 


1 1 


MODEL  FROM  THE  SUFFICIENT  STATISTICS  OF  THE  DATA 

AS  DE  VELOPED  Ill  CHAPTER  2  OF  THIS  THESIS .  LI  USES 
POSTERIOR  DISTRIBUTION  AND  L 2  USES  THE  PRIOR 


DISTRIBUTION  OF  THE  PARAMETERS . 


L  1  J 
L  2  J 

r  a  ] 

L4] 

I.  3  J 
UJ 
[  7 .! 

[  8  J 
I  3  ] 
f  10] 

L 1 1 J 


v  ns s [ !  ]  J  v 
H  S  S  ;  A  E  ;  A  A 
N*-pXi  ;  1  J 


i A B  i AC iL 11 i AD ;DD ; N 


Dl)+-DET{AQ+A  ) 

1 1 l«-(  (  2  xo  VAR  )  *  (  -//  v  2  )  )  x  (  *  (  -SSi  (  2  x  VA R  )  )  ) 

AA*-(  A0+  .  x  MU  0  )  +  A  +  .  *  B 
AB+-(A  0+  .  xMLJQ  )  +  2  x /i  +  .  */i 

AD<-(  (  &)B  )  +  .  x/{  +  .  xB  )  +  ( iy/i/1  )  +  .  x  ( IN  V  /l  0+/1  )  +  .  x/1/1 
L 1 «-  (  *  -  0  .  3  x  ( A  /;  -  ( 4^1 3  )  +  .  x  ( IN  V  A  (J  +  2  x/i  )  +  .  x  /[  D  )  )  x  L 1 1 
L 1  <-  (  L  1  *  (  DET  ( A  0  +/1  x  2  )  )  *  0  .  3  )  x  DD  *  0  .  3 
A  E<-  (  (  §  B  )  +  .  *  A  +  .  *  B  )  -  (  §AA  )  +  .  *  ( I H  V  A  0  +/.  )  +  .  *AA 
AC<-(  AE+  (  (<HMU  0  )  +  .  x/lo  +  .  xMUO  )  ) 

L2+L 1 1 x ( ( DET  A  0  )  *  0  .  3  )  x ( *  0 . 3  xA C ) * ( DD *  0 . 3  ) 


V 
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p 

o 

p 

p 


PROGRAM  'MG  I' 


‘fill  a  PROGRA;!  IE  BASED  OR  MAROjJARDT  '  J  (  1  9  0  3  ) 


A  LG  OR  I  Til  D1  FOR  PARAMETER  ESTIMATION  FRUrl  n 
HON -LINEAR  MODEL. 


VilGTl U jV 
V  MG  T 

[11  LAM*-  o.oi 
[2]  IT*-  0 
1.3]  !«-(  i/./ )  o  .  =  i,7 

L  4 ]  LBB :  -m  0 

[  5  J  PllI*FlJCR  LM 

[G]  X*-DS!1M  LM 

[7]  G*-(§X)  +  .x(Y-PIII) 

[  B  ]  MS*-  (  $  (  7  -  PHI )  )  +  .  x  (  Y  -  P//I ) 

[  9  J  >!«-(  («U')+.  xY  )+/lQ 
[10] 

r.iii  zp-j 

[12]  TL*-  0 

[  1  3  J  L/1J  :  I L*- 1  L  +  l 

[14]  DHL;  IL  J  *-A  L IL  ;  IL  ]  *  "  0  .  5 

L  1  5  J  -+LA  J  x  i  IL  <  M 
[  1  G  J  /!«-//+.  x/l+  .  xD 

[17]  P1<-(P+.  x(.INV(A+I,Al!xI)  )  +  .  x/7  )+.  x/7 

[18]  D2*-(D+.  x(INV(A  +  (LAHilIU)xI )  )+  .  xp  )+  .  x(; 
[10]  Plx-P  +  Pl 

[20]  3  2*-B+D2 

[21] 

[22  1  P+-P1 

[23]  Pnil<-FUCIl  LM 

1-24]  /*-<-/?  2 

[25]  PV,I2*-FUCU  LM 

[  27  ]  Z?l<-($(  Y-PI1I1  )  )  +  .  x  (  Y-Pi/Il  ) 

[  2  3]  S2*-{§{Y-PUI2)  )  +  .  x  (  Y-PI1T  2  ) 

[  2  9  ]  ->■  TIlREx  i  (  (  |  (  +  /  (  BB  -  B  1  )  )  )  <  0  .  0  0  0  1  ) 

[30]  B*-B  1 

[31]  3 

1  32]  TWO  x  i  ,j  1  >  u  L 
[33]  -LBB 
[  3  4  ]  ONE:  LAH*-LAMiNU 
[35  1  *►  LB  B 

[  3  0  J  TWO  :  LAM-LAMxRlJ 
[37]  -LBB 

[  3  B ]  THRU : -FOUR* \SG<G 1 

3  0]  X*-DSNH  LM 

4  0]  VAR+SHE 

[  4 1  ]  CONTINUED  Oil  HE  XT  PAGE 
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[41]  /!«-(  ($A')+.  */)+/!  0 

[42]  C<-  VAR*  (INV  A  ) 

[43]  +FI VE 

[44]  FOUR :  B<-BB 
[43]  VARIES  ill 

[  4  G  ]  C<-VAR*(INV  AA) 

[47]  FI VE : ' PARAMETERS  :  ’ ; B 

[48]  'VARIANCE  :  '  \VAR 

[49]  'VAR-COVAR:  '  \C 


n  PR  OCR  Ail  'EXPO  1* 


a  • EXPD 1 ’  COMPUTES  THE  EXPECTED  CHANGE  IN 
a  ENTROPY  FOR  A  HYPOTHETICAL  EXPERIMENT . 


hexpd i r  n ] v 

V  EXPD  1  ;//;//;  JT  ;  J  ;  //// 

[  1  ]  /7<-l 

[2]  7/1 ;?  7«-y// 7«-2  pO 

[3]  I«-  0 

[4]  w;/-*- 2 

[5]  I7/X«-2 

[6]  7/lP<-10*  3 

[7]  PA+-PW+-  2 

[8]  LAB  :  I  *-I  +  1 

[9]  -*LAT’x  iPI/> 4 0 

[10]  ~>LAS x  \ P/1  >80 

[11]  Ijvo 

[12]  LAC :  IT<-IT+1 

[13]  /10«-JFFFP  IT 

[14]  y//7[IT>S  FPC/V  IT 

[15]  PHI+YHVIIT] 

[16]  x<-B  DSNM  IT 

[17]  VARV[.IT1<-X+.  x(INViAO)  )+.  x($X)xVAR 

[18]  -+LAC*\IT<UM 

[19]  Y11++/PQXY1IV 

[20]  7i4/?7?«-(+/(  VARVxPQ)  )  +  VAR 

[  21  ]  77T«-(®(  7/1/?/?*  (  7/1,77+7/1/?)  )  ) 

[22]  P«-0.  5  x  (  (+/(/'  0  x  (  VTT+  (  (  (y//7-Y//)*2  )  *  7/4/?/? )  )  )  )  ) 

[23]  R  ;  ’  P/1,  Pi/  =  ’  ;  PA  ;  ’  ’  ;  Pi/ 

[24]  ZJ,l^P/t  +  J//y 

[25]  -*F/1P 

[26]  LAS iPW+PW+INX 

[27]  P/1 -<-0.2 

[28]  +L/1P 

[  2  9]  LA  T  :  ->  i  0 

V 
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a  APL  PROGRAMS  FOR  PETROLEUM  RESERVOIR  ESTIMATION . 

a  PROGRAM  'WDR' 

a  ' WDR ’  ESTIMATES  THE  WATER  ENCROACHMENT  TERM .  THIS 
a  PROGRAM  IS  THE  APL  VERSION  OF  THE  FORTRAN  PROGRAM 
a  OBTAINED  FROM  DR.  BENS TEN. 


VWDR [[]]  V 

V  QT+RD  WDR  TD \BEB  \RR \AEA  \F2 \AA \AB \F 1 \F XI \FX2 \LRD 

[I]  -+LAB10*  \RD<33 

[2  J  QT+C  0.2  51233  )  +  (  1  .  503974x21D*0.5)  +  (012  9  55831x?Z)) 

[3]  Q  T+-Q  21+(-0.00  0142643x  TD  *2)  +  ((1.9512£'  7)*TD*3) 

[4]  QT+-QT+(  (“1.1496£,"10)x21I?*4) 

c  5  J  -+LAB2 

[6]  LAB  10:-n0 

[7]  LRD*<*RD 

[8]  FXl<-(  (  (  0.21060396  )  -  0 . 0  1 3  8  30  39  6  xLRD  )  *LRD  ) 

[9]  F*2«-(  (  (  0. 19175244  )-  0 . 0  18  2  0  8  5  0  8  x  LRD  )  x  LRD  ) 

[10]  FX1  +  1 .42555697  +  C ( ( -  1  .  9 6  1 9  3 1 9  ) +FX1 ) *LRD ) 

[II]  FX2*-2 .4  9  81863+(((-l.  7329829)  +  FX  2  )*LRD) 

[12]  F1+(RD*  (  0 . 24  7  786  32  +  (8  .  14  9  72  08£,~6  )  xRD  ))-0. 091510171 

[13]  F2  +  (RD*( ( -1.4646  512£“5  )  +  (  9  .  7  2 0138 3£_ 8  )xRD)  )  +  0  .  0010725209 
[  1  4  J  F1+-  (  FI  *RD  )  -  0 . 1  5  3  5  37  76 

[15]  F  2<-((  (  F2  xRD  )  +  0  . 042167467  )  *RD  )-0. 04562844 

[16]  RR+URD*  2)-1)t2 

[17]  AA+-( *FX 1 )*  2 

[10]  AB*-(  *FX2  )  *  2 

[  19  J  -+LAB  lx  \TD>  0 
[  2  0  J  QT+- 0 

[21 J  +LAB2 
[ 2  2  J  LAB  1 :AATD+AA *TD 

[23]  -+LAB3x  \  AATD  >  112 

[24]  AEA+F1H*AATD) 

[25]  -+LAB  4 

[26]  LAB3-.AEA+-0 

[ 2  7  J  LABA :ABTD+AB*TD 
[28]  -+LAB5*  \ABTD>112 

[29 J  BEB+F2*(*ABTD) 

[30]  6 

[31]  LABb  : BEB-*-0 

[32]  L/ltf6  -.QT+RR-2*  (AEA+BEB  ) 

[33]  LAB2:+\0 

V 
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pi  PROGRAM  'J>RM ' 

pi  'PRM'  ESTIMATES  BOTH  THE  LINEAR  AND  NON- 
pi  LINEAR  PARAMETERS  FROM  THE  MATERIAL 
pi  BALANCE  MODEL  AND  THE  PRODUCTION  HISTORY. 


VPRMl □  ]  V 

V  PRM 

[lj  N+pFET 

[2 J  \N)o  .  =  i 3 

[3]  *[;lJ«-//pl 

[4]  SEARCH 

[  5  ]  XX+  (  i  //  )  o  .  =  1 2 
[6]  XXlilJ^XL;l] 

[  7  J  XXL  ;  2  ]  +-XI  ;  2  j 
L  8  J  C+INV((\XX)+.xXX) 

[9]  BL+C+. x(^XX)t. *FET 
[10J  YC+XX+.xBL 

[ 1 1  J  S++/ (FET-YC)* 2 

[12]  'LINEAR  PARAMETER' \BL 

[13]  ' NONLINEAR  PARAMETERS' \B 

[14]  ^7?i-5*(//-(2+p^[l;])) 

[15] 

[16]  ' LINEAR  VAR-COVAR' 

[17]  C 

[10]  'ERROR  VAR:  '  ; 

V 


pi  PROGRAM  'SEARCH ' 

fl  ' SEARCH '  PERFORMS  A  NON-LINEAR  SEARCH  TO 
pi  ESTIMATE  THE  NON-LINEAR  PARAMETERS . 


VSEARCHlUlV 
V  SEARCH 


[1] 

LAC  :S+-+/(FET-F0FX 

B)*2 

m 

S;' AT  '  \  B ; '  AND 
BB+B 

BL :  '  ;BL 

[4] 

B+-  J 

[  5  J 

-+LAT*  ifl[  1  J  =  0 

[6] 

-+LAC 

[7] 

LAT:B+BB 

V 
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ft  PROGRAM  'FOFX ' 

ft  THIS  PROGRAM  COMPUTES  ALL  THE  TERMS  IN  THE 

ft  MATERIAL  BALANCE  EQUATION  TO  OBTAIN  THE 

ft  DESIGN  MATRIX . 


VFOmCUV 
V  YC+FOFX  B;NN 
[1]  TR+T*Bt  2] 

C  2  J 

[3]  /W<-pT 

[ 4  j  QT+NNpO 

[5]  IL+ 0 

C  6  J  L/M  :IZ>IL+1 

[  7  J  f/Z?/?  ri?[IL] 

[8]  -+LAAx\IL<NN 
[  9  J  J+O 
C  10  J  LAT:J+J+ 1 
CllJ  TERM+O 
[ 1 2  J  I«-0 
[  1  3 J  L/1£:I«-I"+1 

[  1 4  J  2,OA/^2,^/?A/+i7P[«7+5-JJxQTC  JJ 

C  15  j  i/<t7  +  4 

[  1 6  J  XLJ  ;2l<-TERM 
[  1  7  J  +LATx\J<N 

[18]  XZi2l+Xli2l*ET 

[19]  XL  ;  3]«-*[  ;2]*2 
[  2  o  j  /i;r«-(:KJO  +  .x* 

[21]  ZM«-(  \3)°.=i3 

[22]  1+ 0 

[23]  £Itf:I<-I+l 

[24]  Z?A[  J;lJ^/l^[I;r]*'0. 5 
[2  5]  -+SIN*  i  J  <  3 

[26]  i4X«-ZM  +  .  x/lX  +  .  x£>4 

[27]  C+ZM  +  .  x(J/yy  4*)  +  .x£M 

[28]  BL+C+.x(\X)+.*FET 

[29]  YC+X+.*BL 


V 


- 


